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Preface 


The  Intent  of  this  report  is  to  examine  the  robustness  character¬ 
istics  of  a  newly  developed  digital  control  law.  This  study  will 
hopefully  give  insight  to  the  use  of  the  Output  Predictive  Dead-Beat 
Control  law. 
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motivation  and  suggestions  helped  make  this  report  as  complete  as 
possible. 
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Abstract 


A  new  discrete  control  law  (£e£~37"is  implemented  and  examined. 


Robustness  to  noise  and  model  errors  of  the  control  law  as  sample 


rate  varies  is  analyzed.  This  analysis  is  conducted  while  controlling 
several  different  very  lightly  damped,  single  input/single  output  systems 
which  are  representative  of  the  flutter  dynamics  of  the  B-52E  wing  * 


The  control  law  performance  is  different  at  separate  sample  rates. 
A  distinct  range  of  sample  rates  are  found  to  have  a  better  response 
to  noise  than  other  sample  rates.  Another  range  is  found  to  be  more 
robust  when  there  exists  an  error  in  the  models  used  to  calculate  the 


closed  loop  control  law.  When  these  ranges  of  sample  rates  intersect, 
the  robustness  characteristics  at  those  sample  rates  is  found  to  be 


good  with  respect  to  both  noise  and  model  mismatch. 

As  theoretically  predicted^ -<RefT3T"a^relationship  between 
condition  number  of  the  system  Hankel  matrix  and  robustness  seems  to 
exist.  Hence,  these  simulated  results  appear  to  validate  the  theoretical 
results  on  robustness  predicted  by  Reid,  but  on  the  other  hand,  these 
simulated  results  indicate  that  the  total  analysis  of  ^robustness^  is 
a  very  complex  issue  and  cannot,  at  this  point,  be  totally  predicted 
by  such  a  parameter  as  simple  as  the  condition  number  of  the  Hankel 


matrix 


ROBUSTNESS  STUDIES  OF  OUTPUT  PREDICTIVE 


DEAD-BEAT  CONTROL  FOR  WING  FLUTTER  CONTROL  APPLICATIONS 


I  Introduction 


Background 

The  dead-beat  control  law  for  discrete  systems  is  now  well  known 
and  understood.  Simply  stated,  the  dead-beat  control  law  assigns  the 
discrete  time  closed  loop  eigenvalues  to  the  origin.  The  states  will 
be  brought  precisely  to  rest  (assuming  no  additional  disturbances)  in 
no  more  than  n  (system  order)  discrete  steps.  Thus  the  dead-beat  control 
law  has  been  treated  as  an  eigenvalue/eigenvector  assignment  problem 
(Ref  2,  4).  The  dead-beat  control  law  anticipates  the  system  response 
by  feeding  back  all  of  the  system  states. 

Another  approach  for  anticipation  of  the  systems  response  is  to 
actually  predict  the  system  output  into  the  future.  Then,  using  this 
predicted  output  determine  a  control  action  which  forces  the  predicted 
output  identically  to  rest  and  remain  at  rest  with  no  further  required 
input.  Output  Predictive  Dead-Beat  Control  (OPDEC)  uses  this  approach 
in  the  formulation  of  its  control  law. 

In  the  formulation  of  OPDEC,  there  is  no  restriction  on  the 
selection  of  the  discrete  sample  rate.  Theoretical  analysis  of  OPDEC 
produced  an  approach  for  selection  of  an  optimal  sample  rate,  which 
might  enhance  the  systems  over-all  robustness  (Ref  4) . 

Objectives 

This  report  is  concerned  with  the  verification  of  the  robustness 
properties,  with  respect  to  noise  and  model  errors,  of  OPDEC,  and  not 
with  the  underlying  theory  used  in  the  formulation  of  the  control  law. 
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The  initial  objectives  of  this  report  is  to  develop  two  programs  for 
robustness  verification.  The  first  program  finds  the  optimum  sample 
rate  of  a  given  system.  The  second  program  implements  OPDEC  in  a 
simulated  closed  loop  environment  with  noise  and  model  mismatch. 

After  development  of  the  programs  and  selection  of  a  system  to  control, 
verification  of  OPDEC  is  then  performed.  Such  items  as  analysis  of 
robustness  properties  versus’  sample  rate  is  a  major  objective.  Checking 
the  robustness  properties  of  OPDEC  when  using  a  nonminimum  phase  system 
became  another  interesting  area  of  concern.  These  objectives  stated 
above  are  the  prime  areas  of  Investigation  of  this  report. 

Potential  Applications 

This  thesis  is  to  provide  a  basis  for  possible  future  digital 
flight  control  applications  of  OPDEC.  The  OPDEC  concept  appears  to  be 
a  good  candidate  for  digital  flight  control  applications  because  of  its 
"robustness"  properties.  This  characteristic  is  desirable  because  of 
dramatic  changes  in  the  B-52E's  wing  flutter  dynamics  with  changing 
flight  conditions  (altitude  and  airspeed).  The  mathmatlcal  models  that 
describe  the  wing  flutter  dynamics  of  the  B-52E  (Ref  6)  will  be  shown 
later. 
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II  Theory 


An  Important  class  of  control  problems  Is  the  so  called  "tracking 
problem".  The  closed  loop  "tracking  problem"  roughly  Is  the  following. 
For  a  given  reference  variable,  find  an  Input  such  that  the  controlled 
system  output  follows  or  tracks  the  reference.  A  class  of  tracking 
problems  consist  of  those  where  the  reference  variable  Is  a  constant. 
Such  a  problem  is  called  the  regulator  problem  (Ref  1:  Ch  2). 

With  rapid  development  and  miniaturization  of  digital  computers, 
their  use  in  control  systems  has  become  very  common.  A  discrete  con¬ 
troller  that  solves  the  regulator  problem  is  the  "dead-beat"  controller. 
The  "dead-beat"  controller  drives  any  initial  state  to  zero  in  (at  most) 
n  steps,  where  n  is  the  system  order.  The  states,  however,  will  not  be 
driven  to  zero  if  the  output  is  driven  to  zero  (Ref  4 :  Ch  13) .  The 
output  might  have  an  unacceptable  response  between  sample  periods. 

The  Output  Predictive  "Dead-Beat"  Controller  (OPDEC)  derived  by 
Reid  (Ref  3)  has  an  acceptable  response  between  sample  rates  when  the 
output  is  driven  to  zero.  The  OPDEC  control  law  will  drive  the  output 
from  any  initial  point  to  zero  in  at  most  n  discrete  steps.  Because  of 
the  formulation,  the  output  will  also  remain  at  zero  unless  disturbed. 
This  means  the  "states"  of  the  system  are  actually  driven  to  zero  in 
these  n  discrete  steps.  The  needed  sequence  of  control  are  obtained 
from  (Ref  3) 

U(k)  -  -(0.,0.,  ,1) •H^1*Y(k+n/k-l)  (1) 

k  ■  0,1,2,. .. ,n-l 

where  the  output  prediction  vectors 
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[y(k+n/k-l) 

y(k+n+l/k-l) 


Y(k+n/k-l) 


y(k+2n-l/k-l) 


elements  is  the  output  of  the  system  at  discrete  times  k+n,  k+n+1,  . 
k+2n-l  due  to  the  system  states  at  time  k,  or  inputs  up  to  time  k-1, 
and  Hq  is  the  discrete  time  Hankel  matrix  of  size  n.  The  Hankel  matrix, 
shown  below 

h(l)  h(2) . h(n) 

h(2)  h(3) . h(n+l) 

u  _  h(3)  h(4) . h(n+2)  ... 


|_h(n)  h(n+l). . .  h(2n-l}J 

has  only  2n-l  separate  elements.  This  is  helpful  when  trying  to  implement 
OPDEC  on  a  small  computer.  The  components  of  Hq  are  the  discrete  impulse 
response  of  the  system  to  be  controlled. 

According  to  Reid  (Ref  3)  given  a  S1S0  system 

x(t)  -  A*x(t)  +  B*u(t)  (4) 

x(0)  -  Xq  (5) 

with  sampled  output  (sample  time  T) 

y(kT)  -  C*x(kT)  (6) 

that  is  completely  observable  and  completely  controllable,  with  a 
discrete  time  model  of  the  system  (4)-(6)  denoted  as 


x(k+l) 

-  F.x(k)  +  G*u(k) 

(7) 

x(0) 

■ 

(8) 

y(k) 

-  C*x(k) 

(9) 

where 
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G  -  (/0TeAt.dT).B 


(10) 

(ID 


output  can  be  driven  to  zero  In  n  steps  using  OPDEC  to  control  the 
system. 

The  elements  of  the  output  prediction  vector ,  Y(k+n/k-l),  can  be 
found  by 

* 

y(m/k-D  “  C-F®"k-x(k)  (12)  ' 

or  by  the  discrete  conclution  summation 

00 

y(m/k-l)  -  i£oh(m-k+i).u(k-i)  (13) 

Since  equation  (12)  can  be  used  for  prediction  whether  the  open  loop 
system  Is  stable  or  unstable,  this  thesis  will  use  equation  (11)  for 
output  prediction.  This  will  enhance  future  studies  of  OPDEC  to 
controlling  unstable  systems. 

Everything  that  has  been  discussed  so  far  has  not  severely  limited 
us  in  selection  of  the  sample  rate.  The  sample  rate  that  minimizes  the 
condition  number  of  the  Hankel  matrix  yields  controls  with  good  magnitude 
properties  (Ref  4) .  This  thesis  is  investigating  if  this  same  sample 
rate  will  yield  good  robustness  characteristics  as  compared  to  other 
sample  rates.  This  sample  rate,  in  this  report,  is  called  the  "optimal" 
sample  rate. 
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Ill  Investigation 


The  Investigation  into  the  characteristics  of  OPDEC  required  two 
FORTRAN  programs  to  implement  the  algorithm  (See  Appendices  A  and  B). 

The  first  program  implements  a  plot  of  the  condition  number  versus  the 
sample  rate.  This  plot  is  used  to  find  the  optimal  sample  rate  for  the 
sampled  data  controller.'  This  theoretical  "optimum"  occurs  when  the 
reciprocal  condition  number  is  a  maximum  (Ref  3).  To  verify  the  exlstance 
of  an  optimal  sample  rate  with  respect  to  robustness  characteristics, 
the  optimal  sample  rate  and  selected  other  sample  rates  were  then  used 
in  the  implementation  of  OPDEC. 

The  second  program  implements  OPDEC  and  all  the  possible  options 
needed  to  Investigate  the  actual  closed  loop  robustness  properties  of 
OPDEC  (See  Fig  1).  This  program  is  split  into  four  general  sections. 

The  part  initializes  the  program  by  reading  in  the  inputs  and  generation 
of  the  initial  state  vector,  x(o).  Some  of  the  inputs  that  are  read  in 
are  the  true  and  perturbed  state  matrix  equations,  sample  rate  and  logic 
switches.  The  second  part  takes  the  state  vector,  jc(k),  and  predicts 
the  output  in  the  future.  These  predictions  are  put  in  a  vector  format, 
y.  This  is  done  using  either  the  true  model  as  perturbed  model  according 
to  the  switch  one  (SW1)  logic  value.  The  third  part  takes  the  predicted 
output  and  uses  the  Hankel  matrix,  which  it  generates  (Ref  3),  to  find 
the  input,  U.  The  Hankel  matrix  can  be  created  using  the  true  model  or 
a  perturbed  model  according  to  switch  two's  (SW2)  logic  value.  The 
fourth  part  takes  the  input,  U,  and  state  vector  x(k),  and  updates  the 
state  vector  one  sample  time,  x(k+l) .  Also,  at  this  point,  one  can  add 
noise  to  the  input  before  it  is  applied  to  the  true  system.  Also  one  can 
add  noise  directly  to  the  updated  states,  x(ktl),  before  they  are  used  in 
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USING  TRUE 
MODEL  FOR 
PREDICTION 


SING  FALSE 
MODEL  FOR 
REDICTION 


the  next  prediction.  The  noise  added  in  each  case  is  a  zero  Bean  white 
guassian  noise  with  selectable  average  strength.  The  programs  are 
discussed  in  detail  in  appendix  A  and  B. 

Models 

This  thesis  used  basically  two  math  models.  These  models  will  be 
aodified  to  help  answer  some  basic  robustness  questions.  The  models 
are  both  reduced  order  models  of  the  B-52E  wing  flutter  modes  (Ref  6). 

They  both  are  very  lightly  damped.  One  is  a  fourth  order  system  (Table  1) 
and  the  other  is  a  tenth  order  system  (Table  2).  Most  of  the  analysis 
was  done  on  the  fourth  order  system  to  save  computer  simulation  time. 

The  tenth  order  system  was  mainly  used  to  see  effects  of  reduced  order 
models  controlling  larger  systems . 

TABLE  1 

4th  Order  SISO  System 
Open- loop  transfer  function 
8611.7698 

S4+l. 6S3+274 . 075S2+2  79 . 096S+8611. 7698  j 

i 

l 

EIGENVALUES  j 

-.55  +  j  6.0  i 

-.55  -  j  6.0  i 

-.25  +  j  15.4 

-.25  -  j  15.4  j 
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TABLE  2 


10th  Order  SISO  System 
Open-loop  transfer  function 
NUMERATOR 


POLYNOMIAL  ZEROS 


( 

-6607.  )S**9 

( 

.8112E-04)  +  j( 

.19443-01) 

( 

. 2970E+06) S**8 

( 

.8112E-04)  +  j( 

-.1944E-01) 

( 

-.2667E+07)S**7 

< 

-1.843 

)  +  j< 

4.255 

) 

( 

.1405E+09)S**6 

( 

-1.843 

)  +  j( 

-4.255 

) 

( 

-.8689E+08)S**5 

( 

3.343 

)  +  j( 

-15.59 

) 

( 

.2036E+11)S**4 

( 

3.343 

)  +  J( 

15.59 

) 

( 

,5901E+11)S**3 

( 

-2.344 

)  +  JC 

15.48 

) 

( 

.4132E+12)S**2 

( 

-2.344 

)  +  K 

-15.48 

) 

( 

-.4475E+08)S**1 

( 

46.64 

)  +  j( 

0. 

) 

( 

.1561E+09) 

DENOMINATOR 

POLYNOMIAL 

POLES 

( 

1.000  )S**10 

( 

-.9234E-03)  +  j( 

.6579E-01) 

( 

6.626  )S**9 

< 

-.9234E-03)  +  j( 

-.6579E- 

-01) 

( 

948.6  )S**8 

( 

-.6772 

)  +  J( 

-1.053 

) 

( 

4540.  )S**7 

( 

-.6772 

)  +  J( 

1.053 

) 

( 

. 2910E+06) S**6 

( 

-.5446E-01)  +  j( 

15.51 

) 

( 

.8712E+06) S**5 

( 

- . 5446E-01)  +  j( 

-15.51 

) 

( 

. 2919E+08) S**4 

( 

-1.917 

)  +  J( 

16.84 

) 

( 

. 3885E+08) S**3 

( 

-1.917 

)  +  J< 

-16.84 

) 

( 

. 4425E+08) S**2 

( 

-.6631 

)  +  j( 

20.15 

) 

( 

. 2493E+06)S**1 

( 

-.6631 

)  +  i( 

-20.15 

) 

( 

. 1907E+06) 

Questions  to  be  Addressed 

The  emphasis  of  this  thesis  is  to  investigate  robustness  properties 
of  OPDEC  with  respect  to  model  errors  and  noise  addition.  The  obvious 
question  to  ask  is  hov  robust  is  OPDEC,  and  what  factors  enhance  the 
robustness  properties  of  OPDEC.  The  idea  of  an  optimal  sample  rate 
naturally  raises  the  question  of  the  existence  of  an  optimum.  Since 
the  Henkel  matrix  is  instrumental  in  the  selection  of  the  optimum  sample 
rate  and  in  the  usage  of  OPDEC,  it  would  be  interesting  to  see  what 
sensitivities  the  Hankel  matrix  has.  It  would  also  be  interesting  to 
(  examine  the  Issue  of  performance  and  robustness  for  non-minimum  phsse 
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systems.  Controlling  a  system,  with  zeros  In  the  right  half  plane 
In  the  past  has  been  difficult. 


Procedure  and  Results  with  4th  Order  SISO  System 

The  following  results  use  the  same  basic  fourth  order  math  model 
described  previously.  In  some  sections  there  has  been  modifications  to 
this  basic  system.  These  modifications  were  needed  to  answer  some 
specific  questions  and  will  be  discussed  in  each  section. 

The  system  is  assumed  to  start  at  a  random  initial  state  of 


x(o) 


'7.049' 

8.095 

5.007 

9.745^ 


in  phase  variable  coordinates.  For  convenience,  this  coordinate  system 
was  selected  and  due  to  time  considerations,  this  initial  state  was  to 
only  initial  starting  point  used.  The  open  loop  response  of  the  system 
from  this  initial  state  is  shown  in  Figure  2. 


Selection  and  Implementation  of  Sample  Rates 

The  basic  fourth  order  system  was  analyzed  with  the  first  program 
(Appendix  A).  From  this  analysis  of  the  reciprocal  condition  number 
(See  Fig  3)  the  optimum  sample  rate  occurs  at  t  ■  .152,  with  a  maximum 
value  of  1/K  -  .32356. 

Besides  this  optimum  sample  rate,  it  was  desired  to  compare  perfor¬ 
mance  at  several  other  selected  sample  rates.  An  interesting  sample 
rate  to  look  at  is  t  ■  .231.  This  sample  rate  occurs  at  the  relative 
of  the  second  lobe.  This  sample  rate  was  chosen  to  see  if 
relative  optimisation  occurs.  Two  other  sample  rates  chosen  are  .113 
and  .182.  They  were  chosen  because  their  condition  numbers  were  the 
same  as  the  sample  rate  of  .231.  This  was  dona  to  see  if  the  value  of 
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the  condition,  no  natter  where  it  occurs,,  can  be  used  to  predict  the 
robustness  of  the  system.  Another  sample  rate  chosen  was  .954 
because  it  occurs  at  the  relative  maximum  of  the  last  lobe.  It  was 
chosen  because  its  condition  number  is  close  to  the  value  of  the  previous 
three.  The  sample  rates  chosen,  their  reciprocal  condition  number  and 
their  condition  numbers  are  listed  in  Table  3.  The  next  three  sample 
rates  were  chosen  because  of  their  small  condition  numbers.  The  sample 
rates  are  .085,  .214  and  .5.  The  last  three  sample  rates  were  chosen 
because  of  their  very  small  reciprocal  condition  numbers.  The  sample 
rates  are  .04,  .623  and  .835.  These  choices  of  sample  rates  appear  to 
cover  the  broad  range  of  the  condition  number  and  should  provide  enough 
information  for  this  study. 


TABLE  3 

SAMPLE  RATES 
CHOSEN 

1/K 

K 

.152 

.32356 

3.090 

.231 

.10646 

9.392 

.182 

.10714 

9.333 

.113 

.10664 

9.377 

.954 

.09551 

10.469 

.085 

.00925 

108.101 

.214 

.00930 

107.511 

.5 

.00931 

107.366 

.04 

.0004160 

2403.557 

.623 

.00003552 

28151.568 

.835 

.000423 

2363.574 

The  plots  of  the  output  predictive  response  and  the  sequence  of 
control  inputs  for  selected  sample  rates  are  shown  in  Figures  4  thru  16. 
These  plots  were  generated  using  the  true  model  in  both  prediction  and 
control  phases  end  no  noise  added.  Figures  4  and  5  are  the  output 
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response  and  control  Inputs  at  the  optimum  sample  rate  (.152).  The 
output  response  and  control  inputs  for  sample  rates  of  .113,  .182  and 
.231  are  shown  In  Figures  6  and  7.  The  output  response  for  sample  rates 
of  .085,  .214  and  .5  are  shown  in  Figure  8.  The  control  inputs  for 
sample  rates  of  .5  and  .214  are  shown  in  Figure  9.  Due  to  the  large 
scale  needed,  the  control  inputs  for  the  sample  rate  of  .085  is  shown 
separate  in  Figure  10.  The  output  response  and  control  Inputs  for  the 
sample  rate  of  .04  are  shown  in  Figures  11  and  12.  Looking  at  Figures 
10  and  12,  the  control  inputs  for  two  "poorly  conditioned"  sample  rates, 
one  notices  the  inputs  getting  larger  as  the  reciprocal  condition  number 
is  getting  smaller.  This  is  what  was  expected  from  the  theory  (Ref  3). 
Figures  13  and  14  are  the  output  response  and  control  inputs  for  the 
sample  rates  of  .623  and  .954.  The  output  response  and  control  inputs 
for  the  sample  rate  of  .835  is  shown  in  Figures  15  and  16. 

Studying  the  Figures,  one  notices  some  trends.  First,  as  the  sample 
time  increases  the  output  response  starts  fluctuating  more.  This  is 
because  the  sample  rate  is  slower  than  the  natural  frequency  of  the 
system  and  the  system  has  more  time  to  fluctuate.  As  the  sample  rate 
gets  small,  the  value  of  the  control  inputs  increase.  This  is  shown  by 
Figure  7.  Thirdly,  the  condition  number  has  some  effect  on  the  size  of 
the  control  Inputs.  When  both  small  reciprocal  condition  numbers  and 
small  sample  rates  combine,  the  control  input  becomes  huge  (Fig  12). 
Notice  the  magnitude  difference  in  control  input  required  between  sample 
rate  of  .152  and  .04.  This  is  caused  by  the  small  time  increment  the 
system  has  to  achieve  the  dead-beat  response.  For  the  4th  order  system, 
the  system  comes  to  rest  in  4  steps  regardless  of  the  sample  rate.  Thus 
with  a  high  sample  rate  the  system  has  to  work  much  harder  to  drive  the 
states  to  sero  in  the  four  steps.  This  is  the  major  cause  of  magnitude 
difference. 
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SAMPLED  OUTPUT  WITH  NO  NOISE  ADDED 
AT  A  SAMPLE  RATE  OF  .152  SECONDS 
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Fig.  16 


Effect  of  Noise 


This  section  implements  OPDEC  with  input  and/or  state  noise  added. 
This  is  done  for  all  the  sample  rates  selected.  The  average  strengths 
of  the  noises  added  in  either  case  (input  noise  or  state  noise)  are 
.5,  1.0,  1.5  and  2.0.  The  sampled  output  responses,  for  all  selected 
sample  rates,  fell  into  3  categories.  They  are  stable,  conditionally 
stable  and  unstable.  Conditionally  stable  means  the  output  is  varying 
around  the  final  value  and  the  output  is  never  larger  than  it  would  be 
when  no  noise  is  added.  Each  section  has  a  table  that  consolidates  the 
output  result. 

Adding  Input  Noise 

In  this  section,  noise  is  added  eleven  times  to  the  input  while  it 
is  being  input  into  the  true  system.  What  one  would  expect  is  that  as 

i 

the  average  value  of  the  input  noise  is  increased,  the  output  response 

would  become  worse.  This  can  be  seen  in  Fig  17  and  18  which  have 

I 

three  plots  per  figure  at  the  same  sample  rate  with  different  average 
,  input  noises.  Table  4  shows  the  effects  of  input  noise  on  the  sampled 

i 

|  output  at  all  the  selected  sample  rates. 
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SAMPLED  OUTPUT  WITH  INPUT  NOISE  ADDEO  AT  A 
SAMPLE  RATE  OP  .113  SECONDS 
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Pig.  18 
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TABLE  4 


Output  Responses  With  Input  Noise  Added 


SAMPLE  TIME 

OF  CONTROL 

SAMPLE  TIME  OF 
INPUT  NOISE 

OUTPUT  RESPONSE  WITH 
NOISE  STRENGTH  OF 

INPUT 

1/K 

.5 

1.0 

1.5 

2.0 

.4 

.0036 

S 

s 

s 

s 

.000416 

.085 

.0077 

S 

s 

s 

s 

.00925 

.113 

.0102 

s 

s 

cs 

cs 

.10664 

.152 

.0138 

cs 

cs 

cs 

cs 

.32356 

.182 

.0165 

CS 

cs 

cs 

cs 

.10714 

.214 

.0194 

cs 

cs 

cs 

cs 

.0093 

.231 

.021 

cs 

cs 

cs 

cs 

.10646 

.5 

.0454 

cs 

cs 

us 

US 

.00931 

.623 

.0566 

US 

US 

US 

US 

.0000355 

.835 

.0759 

cs 

cs 

US 

US 

.000423 

.954 

.0867 

cs 

us 

US 

US 

.09551 

S  -  STABLE 

CS  -  CONDITIONALLY  STABLE 
US  -  UNSTABLE 

From  Table  4  and  the  control  input  plots  in  the  previous  section, 
one  notices  a  trend.  It  seems  the  effects  of  input  noise  depends  upon 
the  magnitude  of  the  control  inputs.  This  can  be  seen  in  Fig  19  which 
shows  the  sampled  output  with  input  noise  added  at  sample  rates  of  .085 
and  .182.  The  output  response  at  the  "poorer  conditioned"  sample  rate 
(.085)  is  better  than  the  response  at  the  sample  rate  of  .182.  This  is 
because  of  the  magnitude  difference  between  the  respective  control  inputs. 
It  seems  the  larger  the  control  the  less  sensitive  the  output  response 
is  to  input  noise.  It  would  then  be  expected  that  the  output  response 
at  the  sample  rate  of  .835,  would  be  Insensitive  to  input  noise  because 
of  its  large  inputs  relative  to  the  optimal  sample  rate  (Fig  16) . 

Fig  20  is  the  output  response  with  input  noise  added  at  a  sample  rate 
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of  .835.  This  shows  Chat  the  effects  of  input  noise  on  the  output  is 
not  solely  determined  by  the  relative  magnitude  of  the  control  inputs 
or  the  condition  number  before  the  addition  of  the  noise.  It  shows 
that  the  sample  rate  also  effects  the  robustness  of  a  system  that  has 
input  noise  added.  Fig  21  is  an  example  of  the  control  inputs  with 
input  noise  added,  at  a  sample  rate  of  .113  seconds. 

Adding  State  Noise 

In  this  section,  noise  was  added  to  all  the  states  after  they  had 
been  updated  one  sample  step.  This  is  a  simplified  simulation  of  the 
influence  of  incorrect  state  estimation  from  a  Kalman  filter  in  the  closed 
loop  controller.  The  addition  of  noise  occurs  just  before  the  next 
output  prediction.  What  one  would  expect  is  when  the  average  strength 
of  the  noise  is  increased,  the  output  response  becomes  worse.  Fig  22, 

23  and  24  each  have  three  plots  per  figure  at  the  same  sample  rate  with 
different  average  noises  corrupting  the  states.  These  figures  verify 
our  expectations.  Table  5  shows  the  effects  of  state  noise  on  the 
sampled  output  at  all  the  selected  sample  rates. 
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SflHPLEO  OUTPUT  WITH  STATE  NOISE  fiOOEO  AT  A 
SAMPLE  RATE  OF  .113  SECONDS 


36 


Pig.  22 


SAMPLED  OUTPUT  H1TH  STATE  NOISE  ADDED  AT  A 
SAMPLE  RATE  OF  .152  SECONDS 
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Pl«.23 


SAMPLED  OUTPUT  WITH  STATE  NOISE  RODEO  AT  A 
SAMPLE  RATE  OF  .231  SECONDS 
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Fig.  24 
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TABLE  5 


Output  Responses  With  Noise  Corruption  of  the  States 

CONTROL  SAMPLE  OUTPUT  RESPONSE  WITH  STATE 

TIME _ NOISE  STRENGTHS  OF _ 1/K 


.5 

1.0 

1.5 

2.0 

.04 

CS 

CS 

CS 

CS 

.000416 

.085 

s 

CS 

CS 

CS 

.00925 

.113 

s 

S 

CS 

CS 

.10664 

.152 

CS 

CS 

CS 

CS 

.32356 

.182 

CS 

CS 

CS 

CS 

.10714 

.214 

CS 

CS 

CS 

US 

.0093 

.231 

CS 

CS 

CS 

US 

.10646 

.5 

CS 

CS 

US 

US 

.00931 

.623 

CS 

US 

US 

US 

.0000355 

.835 

CS 

CS 

CS 

US 

.000423 

.954 

CS 

CS 

CS 

CS 

.09551 

S  -  STABLE 

CS  -  CONDITIONALLY  STABLE 
US  -  UNSTABLE 

From  Table  5,  one  notices  that  the  stability  of  the  output,  when 
noise  is  added  to  the  states,  is  a  function  of  sample  rate.  With  faster 
sample  rates,  the  system  is  able  to  correct  the  error  induced  by  the 
state  noise.  The  condition  number  is  also  a  factor  as  seen  in  Fig  25. 
The  output  associated  with  the  "poorer  conditioned"  sample  rate  (.214) 
is  not  as  good  as  the  output  at  the  sample  of  .182.  To  best  combat 
state  noise  corruption,  one  would  like  to  use  OPDEC  at  a  fast  sample 
rate  that  has  a  relatively  good  condition  value.  For  the  model  used, 
a  sample  rate  of  .113  seems  to  have  the  best  results  (Fig  22). 

Input  and  State  Noise  Added 

This  section  concerns  itself  with  the  addition  of  both  input  and 
•  state  noise.  Looking  at  Tables  4  and  5,  and  trying  to  anticipate  which 
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sample  rate  would  give  the  best  output  results  when  both  input  and 
state  noise  is  added,  one  would  probably  choose  the  following  sample 
rates  of  .04,  .085,  .113,  .152,  .182  and  .231.  After  a  multitude  of 
computer  runs,  implementing  OPDEC  with  both  noises  added  and  at  all  the 
selected  sample  rates,  the  expected  sample  rates  show  the  best  robustness. 
Figures  26,  27,  and  28  show  the  sampled  output  with  both  input  and  state 
noise  added,  and  as  one  expects,  .113  is  the  most  robust. 

This  section  has  investigated  the  best  way  to  combat  noise  in  a 
system  when  using  OPDEC  for  control.  The  conclusion  appears  to  be,  to 
run  the  system  at  a  fast  sample  rate  that  is  also  well  conditioned.  The 
system  seems  to  be  able  to  correct  the  induced  errors  caused  by  noise 
when  it  is  at  a  faster  sample  rate.  In  effect  the  errors  can  be 
"negated"  before  they  significantly  degrade  the  system  response  when 
operating  at  the  faster  sample  rate. 

Effects  of  Model  Mismatch 

The  program  that  implements  OPDEC  also  has  the  option  to  use  a 
perturbed  model.  This  perturbed  model  can  be  used  in  either  the  predic¬ 
tion  phase,  control  phase  or  both  phases  simultaneously.  A  mathmatical 
model  can  only  approximate  a  system's  response,  but  this  mathmatical 
model  is  used  to  try  and  control  the  system.  Since  the  model  is  in 
error  the  controller  must  be  able  to  cope  with  this  error  without  causing 
the  system  to  go  unstable.  The  perturbed  model  analysis  in  this  section 
tries  to  simulate  these  model  mismatching  errors.  Using  this  perturbed 
model  analysis  enables  a  person  to  see  which  phase  (prediction  or  control) 
is  the  most  sensitive  to  model  errors. 

To  generate  the  "perturbed  model”  the  basic  fourth  order  model 
eigenvalues  are  shifted  in  the  model  versus  the  "true  system".  The 
reasons  for  perturbing  the  eigenvalues  instead  of  the  Hankel  matrix 
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Fig.  28 
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perturbations  are  not  well  understood.  These  perturbed  eigenvalues 
are  then  used  to  create  the  perturbed  model.  Specifically,  two  perturbed 
models  are  used  in  the  section  to  demonstrate  robustness  characteristics 
of  OPDEC.  The  transfer  function  for  the  two  perturbed  systems  are: 

10%  -perturbation 
+  8474 

S4+l . 66S3+314 . 1S2+360 . 6S+8474 

EIGENVALUES 

-.605  +  j  5.4 
-.605  +  j  5.4 
-.225  +  j  16.94 
-.225  +  j  16.94 

15%  perturbation 
+  8285 

S4+l. 69S3+340 .6S2+408S+8285 

EIGENVALUES 

-.6325  +  j  5.1 
-.6325  -  j  5.1 
-.2125  +  j  17.71 
-.2125  -  j  17.71 

OPDEC' s  program  was  then  implemented  at  all  the  selected  sample 
rates,  using  the  perturbed  model  in  the  following  fashion: 

1.  The  perturbed  model  was  used  in  the  prediction  phase  and  the 
true  model  was  used  in  the  control  phase.  (False  prediction, 
true  control) 

2.  The  true  model  was  used  in  the  prediction  phase  and  the 
perturbed  model  was  used  in  the  control  phase.  (True  prediction, 
false  control) 

3.  The  perturbed  model  was  used  in  both  prediction  and  control 
phase.  (False  prediction,  false  control) 
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The  sampled  output  responses,  for  all  selected  sample  rates,  were  stable. 

According  to  Reid  (Ref  3)  the  norm  of  the  difference  between  the 
true  Hankel  matrix  and  the  error  Hankel  matrix  (E  *  Hj  -  He)  gives  an 
indication  of  the  stability  of  the  output.  Tables  7  and  8  show  the 
selected  sample  rates,  the  condition  number,  the  norm  of  the  true  Hankel 
matrix  { ( Ht I ! >  the  norm  of  the  difference  between  the  true  Hankel  matrix 
and  the  error  Hankel  matrix,  | |e| |,  the  norm  of  the  E  matrix  divided  by 
the  norm  of  the  true  Hankel  matrix,  and  the  condition  number  times  this 
value  for  the  10%  and  15%  error  models  used. 


TABLE  6 

10%  Perturbed  Model 


DELT 

K 

l|Htl| 

INI 

1 1 E | | / | i Ht | | 

K.||E||/||Ht|| 

.040 

2403.555 

.652 

.078 

.120 

289.085 

.085 

108.102 

1.562 

.401 

.257 

27.754 

.113 

9.055 

1.566 

.553 

.353 

3.197 

.152 

3.091 

1.718 

.559 

.325 

1.006 

.182 

9.334 

2.086 

.937 

.449 

4.191 

.214 

107.511 

2.366 

1.121 

.474 

50.957 

.231 

9.393 

2.406 

.874 

.363 

3.412 

.500 

107.367 

4.001 

2.543 

.635 

68.227 

.623 

28151.136 

2.012 

2.941 

1.462 

41155.396 

.835 

2363.574 

1.350 

.911 

.675 

1594.370 

.954 

10.470 

.836 

.682 

.816 

8.543 
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TABLE  7 


15%  Perturbed  Model 


DELT 

K 

1  1  Ht  |  | 

llEll 

1 |E| |/| |Ht| | 

K*  | lE[ I/|  |Ht| | 

.040 

2403.555 

.625 

.123 

.189 

454.633 

.085 

108.102 

1.562 

.566 

.362 

39.147 

.113 

9.055 

1.566 

.750 

.479 

4.336 

.152 

3.091 

1.718 

.750 

.436 

1.349 

.182 

9.334 

2.086 

1.226 

.588 

5.486 

.214 

107.511 

2.366 

1.287 

.544 

58.498 

.231 

9.393 

2.406 

1.288 

.535 

5.026 

.500 

107.367 

4.001 

3.589 

.897 

96.315 

.623 

28151.136 

2.012 

3.141 

1.561 

43943.039 

.835 

2363.574 

1.350 

1.104 

.817 

1932.136 

.954 

10.470 

.836 

.903 

1.081 

11.314 

Looking  at 

the  tables. 

at  the  same  sample  rate 

the  norm  of  E 

increases  as  the 

amount  of  perturbation  increases. 

If  the  norm  of  E 

is 

less  than  one 

,  the  output 

response 

is  suppose  to 

be  stable.  This 

section  will  try  and  verify  this  theoretical  result  by  looking  at  the 
output  response  at  all  of  the  sample  rates  selected  using  the  perturbed 
models  and  seeing  if  the  output  is  stable  or  unstable.  Table  8  shows 
the  effects  which  the  10%  perturbed  model  has  on  the  sampled  output. 
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TABLE  8 


Output  Responses  With  the  10Z  Perturbed  Model 


OUTPUT  RESPONSE  WITH 

TRUE  PREDICTION  FALSE  PREDICTION  FALSE  PREDICTION 


SAMPLE  RATE 

FALSE  CONTROL 

TRUE  CONTROL 

FALSE  C01 

.04 

S 

S 

S 

.085 

S 

S 

S 

.113 

S* 

S* 

S 

.152 

S 

S* 

S 

.182 

S* 

US 

S* 

.214 

US 

US 

US 

.231 

US 

US 

s 

.5 

US 

US 

US 

.623 

US 

US 

US 

.835 

US 

US 

US 

.954 

US 

US 

s 

The  time  needed  to  reach  the  final  value  was  greater  than  20  times 
the  sample  rate  used. 


The  output  responses  using  the  10Z  perturbed  model,  for  sample  rates 
of  .04,  .085,  .113,  and  .152  are  in  Figs  29,  30,  31  and  32  respectively. 
From  these  figures  one  can  see  that  the  prediction  phase  is  more  sensi¬ 
tive  than  the  control  phase  to  model  error.  But  what  is  really  inter¬ 
esting  is  the  output  response  is  better  when  the  error  model  is  used 
in  both  prediction  and  control  phases. 

Table  9  shows  the  effects  the  15Z  perturbed  model  has  on  the 
sampled  output. 


48 
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SAMPLED  OUTPUT  WITH  HO  NOISE  RDOEO  USING  THE 

10*  PERTURBED  MODEL  AT  A  SAMPLE  RATE  OF  .085  SECONDS 
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SAHPLEO  OUTPUT  WITH  NO  NOISE  ADDED  USING  THE 

10  %  PERTURBED  MODEL  AT  A  SAMPLE  RATE  OF  .113  SECONDS 
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Pi«.  31 
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SAMPLED  OUTPUT  WITH  NO  NOISE  flDDEO  USING  THE 

lOfC  PERTURBED  MODEL  AT  A  SAMPLE  RATE  OP  .152  SECONDS 
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Pig.  32 
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TABLE  9 


Output  Responses  With  15%  Perturbed  Model 
OUTPUT  RESPONSES  WITH 


SAMPLE  RATE 

TRUE  PREDICTION 
FALSE  CONTROL 

FALSE  PREDICTION 
TRUE  CONTROL 

FALSE 

FALSE 

PREDICTION 

CONTROL 

.04 

S 

S 

S 

.085 

S 

S 

S 

.113 

US 

US 

S 

.152 

US 

US 

US 

.182 

US 

US 

US 

.214 

US 

US 

US 

.231 

US 

US 

US 

.5 

US 

US 

US 

.623 

US 

US 

US 

.835 

US 

US 

US 

.954 

US 

US 

US 

S  -  STABLE 
US  -  UNSTABLE 

The  models  errors  are  almost  too  great  to  overcome  when  the  perturb¬ 
ation  is  15%.  The  output  responses  using  the  15%  perturbed  model  for 
sample  rates  of  .085  and  .113  are  in  Fig  33  and  34.  These  sample  rates, 
as  predicted  by  Table  7,  are  stable.  It  seems  that  the  norm  of  the 
difference  between  the  two  Hankel  matrix  and  the  error  Hankel  matrix 
if  less  than  one  indicates  stability. 

Combination  of  Noise  and  Model  Errors 

This  section  investigates  the  efforts  of  input  and  state  noise 
when  erroneous  models  are  used  to  control  the  system.  OPDEC's  program 
was  then  implemented  with  the  perturbed  model  used  in  both  prediction 
and  control  phases  with  just  input  noise,  then  just  state  noise,  and 
finally  both  input  and  state  noise.  The  justification  for  implementing 
OPDEC  in  this  fashion  is  that  one  will  not  be  using  different  models  for 


the  production  phase  as  compared  to  the  control  phase.  The  same  model 
will  be  used  in  both  phases. 

OPDEC’ s  program  was  implemented  using  the  10%  perturbed  model, 
input  noise  was  added  first,  then  state  noise,  then  the  combination  of 
the  two  noises.  This  was  done  for  all  the  sample  rates  selected,  except 
the  sample  rates  which  by  Table  4  were  unstable.  Shown  in  Figures  35 
and  36  is  the  sampled  outputs  using  the  10%  perturbed  model  with  input 
noise  added.  The  sampled  outputs  with  state  noise  added  is  shown  in 
Figures  37  and  38.  The  sampled  outputs  with  both  state  and  input  noise 
added  can  be  seen  in  Figures  39  and  40. 

As  expected,  the  faster  sample  rates  which  have  larger  controls 
are  insensitive  to  input  disturbances,  but  the  addition  of  state  noise 
has  a  more  severe  effect  on  the  output.  The  severity  is  great  enough  to 
make  the  output  response  worse  than  the  response  at  a  sample  rate  of  .113. 
When  using  a  perturbed  model,  the  influence  of  state  noise  becomes  more 
severe.  This  is  because  of  OPDEC's  sensitivity  in  the  prediction  phase 
of  its  control  law. 

Hankel  Matrix  Sensitivity  Analysis 

The  use  of  OPDEC  requires  working  with  the  Hankel  matrix.  Since 
the  Hankel  matrix  is  such  a  vital  part  of  OPDEC,  the  analysis  of  the 
Hankel  matrix  under  variations  in  the  system  pole/zero  structure  will 
give  some  indication  of  the  overall  robustness  of  OPDEC.  Not  only  is 
the  Hankel  matrix  an  intrinsic  part  of  control  computation  in  OPDEC, 
it  is  also  an  Important  part  in  the  selection  of  an  "optimal"  sample 
rate.  The  approach  taken  in  this  report  is  to  see  how  sensitive  the 
Hankel  matrix  is  to  pole  and  zero  placement. 

First  the  dominate  pair  of  eigenvalues  in  the  fourth  order  system 
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Fig.  38 
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Fig.  40 
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(-.25  +  j  15.4)  are  perturbed  slightly,  then  combined  with  the  two  other 
original  eigenvalues  to  come  up  with  a  new  perturbed  system.  This  system 
is  then  used  in  the  first  FORTRAN  program  (Appendix  A)  to  see  the  effects 
which  pole  placement  has  upon  the  condition  number.  The  resulting  con¬ 
dition  number  versus  sample  rate  plot  is  then  compared  to  the  original 
system  condition  number  plot.  The  original  system  condition  number  plot 
is  again  shown  in  Figure  41.  Figure  42  and  Table  10  show  the  perturbed 
eigenvalues  used  in  the  sensitivity  analysis.  Table  10  also  has  the 
optimal  sample  rate  and  one  over  condition  number  for  each  perturbed 
eigenvalue . 


TABLE  10 

MAXIMUM  RECIPROCAL 


PERTURBED 

EIGENVALUE 

"OPTIMAL" 
SAMPLE  RATE 

CONDITION  NUMBER  @ 
"OPTIMUM"  SAMPLE  RATE 

Wn 

DAMPING 

RATIO 

-.25  +  j  15.4! 

.152 

.32356 

15.402 

.016 

-.25  +  j  20 

.118 

.23109 

20.001 

.0125 

-.25  +  j  17.4 

.136 

.2881 

17.401 

.014 

-.25  +  j  12.7 

.174 

.2963 

12.702 

.019 

-.25  +  j  10.0$ 

.215 

.21664 

10.003 

.024 

-.05  +  j  15.4 

.151 

.3581 

15.4 

.0006 

-.15  +  j  15.4 

.152 

.3399 

15.4 

.0097 

-.35  +  j  15.4 

.152 

.305 

15.4 

.0227 

-.45  +  j  15.4 

.152 

.2877 

15.403 

.0292 

-.2  +  j  12.32* 

.177 

.2891 

12.321 

.016 

-.3  +  j  18.48* 

.128 

.2585 

18.482 

.016 

!  The  original  system 

*  The  perturbed  system  created  using  these  eigenvalues  has  the 
same  damping  ratio  as  the  original  system 

$  Local  optimal  sample  rate  only;  i.e.,  there  is  another  peak 

Looking  at  Figures  43  thru  52  one  notices  a  trend  appearing.  When 
the  real  part  of  the  eigenvalue  is  perturbed  and  the  imaginary  part  is 
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not,  the  effects  upon  the  optimal  sample  rate  is  minimal.  Conversely, 
when  the  real  part  is  not  perturbed  and  the  imaginary  part  is  perturbed, 
the  shift  of  the  optimal  sample  rate  is  considerable.  This  seems  reason¬ 
able  because  the  frequency  of  a  system  is  very  dependent  upon  the 
imaginary  portion  of  the  eigenvalues  in  a  lightly  damped  system.  This 
analysis  shows  the  Hankel  matrix  condition  number  being  sensitive  to 
the  frequency,  Wn,  of  a  system,  but  insensitive  to  the  damping  ratio  of 
a  system. 

The  next  area  of  concern  is  the  effect  a  system  zero  would  have  on 
the  Hankel  matrix.  A  system  zero  is  added  to  the  basic  fourth  order 
system  giving  us  a  new  system,  shown  below. 

(S+«) 8611. 7698 

S4+1.6S3+274.075S2+279. 096+8611. 7698 

Alpha  is  varied  from  one  to  minus  one  to  see  the  sensitivities  in 
the  Hankel  matrix.  Table  11  shows  the  optimum  sample  rate  and  the  corre¬ 
sponding  reciprocal  condition  value  for  various  values  of  alpha. 

TABLE  11 

"OPTIMAL"  RECIPROCAL 


ALPHA 

SAMPLE  RATE 

CONDITION  NUMBER 

1.0 

.147 

.6234 

.3 

.156 

.6336 

.2 

.156 

.6352 

0.0 

.156 

.6395 

-.2 

.156 

.6438 

-.3 

.156 

.6456 

-1.0 

.153 

.6526 

The  original  system  condition  number  plot  is  again  in  Figure  53. 

For  alpha's  of  .3,  0.0,  and  -.3,  plots  of  the  reciprocal  condition  number 
versus  sample  rate  are  shown  in  Figures  54,  55  and  56,  respectively. 
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The  most  obvious  change  is  the  value  of  the  maximum  reciprocal  condition 
number.  The  value  has  about  doubled  from  .3  for  the  no  zero  system  to  .6 
for  the  system  with  a  zero.  Also  observed  is  that  variation  on  the  optimal 
sample  rate  itself  is  very  small  and  it  changes  very  little  as  the  zero 
shifts,  even  when  the  zero  is  in  the  right  half  plane.  The  system  is 
then  a  nonminimum  phase  system,  but  there  seems  to  be  little  influence 
upon  the  characteristics  of  the  conditioning  of  the  Hankel  matrix.  If 
anything,  the  maximum  reciprocal  condition  number  is  increasing  in 
value  as  we  proceed  further  into  the  right  half  plane!  If  the  condition 
numbers  value  is  any  indication  of  robustness,  this  indicates  that  a 
system  with  a  zero  will  be  more  robust  than  the  same  system  without  a 
zero.  From  this  analysis,  using  a  lightly  damped  system,  the  Hankel 
matrix  is  sensitive  to  the  pole  locations  of  the  system,  but  it  is 
relatively  insensitive  to  the  location  of  the  system  zeros. 

Minimum  Phase  versus  Nonminimum  Phase 

This  section  examines  whether  or  not  a  nonminimum  phase  system 
will  present  any  problems  when  being  controlled  by  OPDEC.  Two  systems, 
one  being  nonminimum  phase,  the  other  minimum  phase,  are  implemented 
under  conditions  of  model  errors  and  input  and/or  state  noise.  The 
two  systems  are: 

System  A  (Minimum  Phase) 

8611.7698(S+.3) 

S4+l. 6S3+274 . 075S2+279. 096S+8611 . 7698 

System  B  (Nonminimum  Phase) 

8611.7698(S-.3) 

S4+1.6S3+274.075S2+279.096S+8611.7698 
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From  Table  11,  both  System  A  and  B  have  an  optimal  sample  rate  of 
.156  seconds.  Their  sampled  outputs  and  control  inputs  at  the  optimal 
sample  is  shown  in  Figures  57  and  58.  Surprisingly,  the  sampled  outputs 
and  control  inputs  are  exactly  the  same  for  a  minimum  phase  system  and 
a  nonminimum  phase  system.  This  is  caused  by  the  zeros  in  both  systems 
having  the  same  magnitude.  Input  noise,  then  state  noise,  and  finally 
input  and  state  noise  are  added  to  the  system.  The  sampled  output  with 
input  noise  added  at  a  sample  rate  of  .156  can  be  seen  in  Figure  59. 

Again  the  output  is  exactly  the  same  for  a  nonminimum  phase  system  as 
for  a  minimum  phase  system.  This  fact  is  present  when  state  noise  is 
added  and  also  when  both  types  of  noises  are  added.  Other  sample  rates 
are  tried  and  the  outputs  are  identical  for  system  A  and  B.  Similar 
to  previous  results,  a  faster  than  optimum  sample  rate  gives  better 
robustness  characteristics  when  noise  is  disrupting  the  system  (Figures 
60  and  61) . 

The  next  area  of  concern  is  if  system  A  and  system  B  have  the  same 
properties  when  an  error  model  is  used  for  control.  The  two  10%  perturbed 
models  are: 

System  A' 

8474(S+.3) 

S4+1.66S3+317.1S2+360.6S+8474 
System  B' 

8474(S-.3) 

S4+1.66S3+317.1S2+360.6S+8474 

The  denominator  of  the  perturbed  systems  is  the  same  as  the  previous  10% 
perturbed  models  denominator  used  earlier  in  this  report.  The  output 
response  and  control  inputs  when  using  the  perturbed  models  to  control 
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the  true  systea  are  identical.  Figure  62  shows  the  sampled  outputs  at 
saaple  rates  of  .113  and  .152  using  this  perturbed  model. 

Froa  this  analysis,  OPDEC  has  no  probleas  controlling  a  nonainiaua 
phase  systea.  Also,  the  robustness  of  OPDEC  is  not  degraded  when  this 
type  systea  is  being  controlled.  OPDEC  is  fairly  insensitive  to  zero 
placement. 

This  insensitivity  to  the  systea  zero  location  was  interesting 
enough  to  try  and  control  systea  B  using  systea  A  and  A'  as  the  perturbed 
model.  When  systea  A  was  used,  the  output  response  looked  exactly  like 
the  original  dead-beat  response  (Figure  57).  When  systea  A'  was  used, 
the  output  response  looked  exactly  like  Figure  62.  The  next  step  was 
to  try  and  use  the  original  fourth  order  SISO  systea  with  no  zeros  to 
control  system  B.  The  systems  response  was  unstable.  What  this  analysis 
shows  is  that  the  zero  location  is  not  a  factor  in  using  OPDEC  as  a 
regulator. 

Discussion  of  4th  Order  Results 

From  this  fourth  order  analysis,  we  see  that  saaple  rates  faster 
than  the  optimum  have  better  output  characteristics  when  noise  is 
corrupting  the  system.  Also,  some  of  the  faster  saaple  rates  seem  to 
have  better  robustness  characteristics  overall.  The  condition  number 
is  a  helpful  guide  in  selecting  general  areas  of  saaple  rates  where 
the  system  is  more  robust.  Using  the  condition  number  to  find  a 
specific  "optimal"  sample  rate,  without  doing  the  closed  loop  robustness 
simulation  studies,  is  definitely  the  wrong  approach.  Robustness  end 
closed  loop  system  performance  is  just  too  complex  of  an  issue  to  hope 
that  there  would  be  a  single  panacea  for  all  aspects  of  the  problem. 

Areas  of  high  sensitivity  of  OPDCC  seem  to  be  the  poles  of  the 
system  being  controlled,  but  omc  appears  to  be  insensitive  to  systems 
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zero  location.  As  far  as  the  location  of  the  poles  are  concerned, 

OPDEC  is  sore  sensitive  to  the  "frequency"  of  the  pole  as  compared  to 
the  damping  ratio. 

Procedure  and  Results  with  10th  Order  SISO  System 

In  this  section,  OPDEC  is  implemented  using  the  10th  order  system 
discussed  earlier.  Figure  63  shows  the  reciprocal  condition  number 
versus  sample  rate.  The  "optimum"  sample  rate  is  .523  and  the  reciprocal 
condition  number  is  0.000158.  The  sampled  output  at  a  rate  of  .523  is 
shown  in  Figure  64.  The  output  is  very  oscillitory  but,  as  the  theory 
states,  in  ten  steps  of  the  dead-beat  control  the  output  and  the  states 
are  brought  identically  to  zero. 

What  was  tried  at  this  point  was  to  control  the  10th  order  system 
with  a  smaller  ordered  system.  The  full  ordered  Hankel  matrix  was 
created  (10  X  10).  Then  an  option  was  exercised  in  the  program  that 
lets  the  user  reduce  the  size  of  the  Hankel  matrix  just  before  it  is 
Inverted  for  use  in  closed  loop  control.  Then  the  reduced  order  Hankel 
matrix  is  Inverted  and  used  to  determine  a  control  input  using  OPDEC 's 
control  law.  When  trying  to  control  the  true  10th  order  system  with  a 
9th  order  system  the  output  for  sample  rate  of  .523  is  shown  in  Figure  65. 
For  the  same  sample  rate  the  order  of  the  Hankel  matrix  was  reduced  to 
an  8th  order  system.  The  sampled  output  can  be  seen  in  Figure  66. 

Some  other  sample  rates  were  chosen  and  the  same  procedure  was  used. 

The  outputs  for  all  the  other  selected  sample  rates  went  unstable, 
even  for  the  9th  order  system  controlling  the  full  10th  order  system. 

The  condition  number  is  helpful  in  determining  a  sample  rate  in  which 
a  reduced  order  model  la  being  used  to  determine  control  inputs. 
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IV  Conclusions  and  Recommendations 

From  this  analysis  one  can  see  that  the  sample  rate  can  enhance 
the  robustness  properties  of  OPDEC.  The  sample  rate  which  minimizes 
the  condition  number  does  yield  controls  with  good  magnitude  properties. 

If  anything,  this  fact  causes  the  response  to  be  more  sensitive  to 
input  noise  if  the  strength  of  the  noise  is  held  constant.  Thus,  the 
so  called  optimum  sample  rate  is  not  necessarily  the  same  sample  rate 
that  gives  the  best  enhancement  to  OPDEC's  robustness. 

The  group  of  sample  rates  which  has  good  robustness  properties 
when  input  noise  was  the  only  disruptive  element  is  not  the  same  group 
of  sample  rates  that  exhibit  acceptable  robustness  when  model  mismatch 
exists.  This  means  that  the  overall  characteristics  of  OPDEC  can  be 
tailored  for  specific  operating  environments,  by  just  selecting  the 
proper  sample  rate.  The  problem  then  becomes  one  of  selection  of 
the  proper  sample  rate  to  "tune"  the  algorithm  for  the  particular  set 
of  noise  environment  and  model  errors  which  is  expected  to  be  encountered. 
There  does  not  appear  to  be  an  analytic  method  for  this  proper  sample  rate. 
Rather,  the  robustness  performance  appears  to  be  an  issue  which  is  best; 
analyzed  through  simulation  studies. 

OPDEC's  sensitivity  to  the  frequency  of  the  system  and  its  opera¬ 
ting  sample  rate  could  be  caused  from  the  selection  of  the  system  to 
control.  It  is  recommended  that  future  studies  see  if  OPDEC  has  the 
same  characteristics  when  implemented  on  other  system,  particularly  on 
those  which  are  not  so  lightly  damped  as  the  one  studied  here.  Another 
area  needing  further  investigation  is  in  the  optimization  of  sample  rate. 
The  flexibility  of  OPDEC  seems  to  make  it  a  viable  control  technique. 
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Appendix  A 


Condition  Number  Program 

This  Appendix  contains  the  computer  listing  and  a  users  guide  for 
the  program  that  finds  the  optimal  sample  rate.  The  optimal  sample  rate 
occurs  when  one  over  the  condition  number  is  a  maximum.  This  program 
uses  many  of  the  same  subroutines  as  the  program  that  implements  OPDEC 
(Appendix  B) .  To  save  space,  these  subroutines  will  be  discussed  in 
Appendix  B. 

The  user  must  supply  the  program  with  the  following  information  on 
data  cards.  All  data  is  read  in  using  an  unformated  read  statement. 

1.  System  size  n,  n  is  an  intiger  value  and  must  be  less 
than  or  equal  to  ten. 

2.  Starting  sample  rate,  DELT. 

3.  Final  sample  rate,  DF. 

4.  A,  B,  and  C  from  state  matrix  equation 

x  **  Ax  +  Bu 
y(t)  "  Cx(t) 

A,  B,  and  C  re  all  n  by  n  matrices  and  must  be  discret- 
izable  over  sample  rate  range  user  supplied.  A  must  also 
be  inver table. 

Below  is  a  brief  outline  of  the  steps  the  program  follows: 

1.  Reads  in  data. 

2.  Sets  DET  -  0.001  and  counter  IS  -  0. 

3.  Increments  starting  sample  rate, 

DELT  -  DELT  +  DET 

4.  Checks  to  see  if  DELT  >  DF;  if  so,  goes  to  step  12. 
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5.  Discretizes  A  and  B  at  sample  rate  DELI. 

6.  Increments  counter  IS  *  IS  +  1. 

7.  Creates  Hankel  matrix. 

8.  Finds  Singular  values  of  Hankel  matrix. 

9.  Find  one  over  condition  number, 

YS(IS)  -  Q(n)/Q(l) 

Q(n)  is  minimum  singular  value, 

Q(l)  is  maximum  singular  value. 

10.  Sets  TSS(IS)  -  DELT. 

11.  Go  to  step  3. 

12.  Search  reciprocal  condition  number  array  (YS)  for 

maximum  reciptrocal  condition  number  and  output  this  number. 

13.  Find  the  sample  rate  at  which  the  maximum  reciprocal  condi¬ 
tion  number  occurs  and  output  this  number. 

14.  Plot  array  YS  vs  array  TSS. 

15.  Stop. 

The  program  uses  two  subroutines  from  the  IMSL  subroutine  package. 

The  two  routines  are: 

1.  LINV2F 

2.  LSVDF 

Description  of  Subroutines 

Subroutine  COND(F2,G2,C,N,IDIM,Q) 

This  subroutine  takes  the  discretized  system  and  creates  the  Hankel 
matrix  from  the  observability  and  controllability  pairs.  Then  it  does  a 
singular  value  decomposition  of  the  Hankel  matrix.  The  singular  values 
are  then  output  in  the  subscripted  array  Q. 
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Inputs 

1.  F2,  G2,  C,  discretized  system. 

2.  n,  system  size. 

3.  IDIU,  Initial  dimensionalization. 

Outputs 

1.  Q,  n  by  1  array  containing  the  ordered  singular  values 

Q(l)  maximum  singular  value, 

Q(n)  minimum  singular  value. 

Subroutines  EFT,  HGRAGH,  and  VGRAGH  are  described  in  Appendix  B. 
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PROGRAM  THESIS (INPUT, OUT PUT* TAPES *INPUT ,TAPE6=0UTPUr ,  PLOT) 
01 ME NS E ON  ATCiOfiOJ ,BT (1 0,1b ) ,CT(1Q, 10) 

DIMENSION  Q  (10) 

DIMENSION  TSS(2500)  ,YS(2500> 

DIMENSION  IYSS (17) 

DIMENSION  FT(1Q,10),GT(10»10) 

OATA  IfSS(l)/20HONE  OVER  CONDITION  / 

OATA  If SS (3)/2QH  NUMBER  VS  SAMPLE  / 

OATA  IfSS(5)/20H  RATE  / 

DATA  If SS (7J/20H  / 

OATA  If  SS  (9)/2QH  SAMPLE  RATE  / 

OATA  IYSS (11)/ 20 HONE  OVER'  CONDITION  #/ 

OATA  If SS(13)/AOHON£  OVER  CONDITION  NUMBER  VS  SAMPLE  RATE/ 
IS=0 

C  READ ING  IN  SIZE  OF  Mf  SYSTEM  N 
C  N  MUST  BE  t.£SS  THAN  OR  EQUAL  TO  10 
READ*,N 
IDIM=13 
IPASs0 

C  READ  IN  INITIAL  OELT 
REAO*,OELT 
OET-O.OOl 
C  READ  FINAL  T 
READ*, OF 

REAO*, ( (AT (I, J) , J=lf N) ,I*1,N) 

REAO*, ((BT(I,J),J=1,N>, 1*1, N) 

REAO*,  (  (CT (I*  J)  , J*1»N)  »I=1»N) 

12  CONTINJE 

DELT*DET*f)ELT 
IF(OELT.GT.OF)  GOTO  96 
IS*IS+l 
IPAS»I3AS*1 

C  GENERATE  Mf  OISCRETE  F  AND  G  MATRIX  FOR  TRUTH  AND  MODEL* 

Ml*55 

CALL  EsrT(ATfBT»N»IOIMfFT  ,GT»Mi,D£LT> 

CALL  COMO(FT,GT,CT,N,IOIM,Q) 

TSS(IS) *OELT 
YS(IS)-3(N)  /Q  ( 1) 

IF(IPAS t GE» 10) IPAS*0 
IFdPAS *NE. 0)  GOTO  12 
PRINT*,"  " 

PRINT*,"  SAMPLE  RATE  IS  ",TSS(IS> 

PRINT*,"  " 

PRINT*,"  ONE  OVER  CONDITION  NUMBER  IS" 

PRINT*,"  1/K*", YS (IS) 

GOTO  12 

96  CONTINJE 
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290 

269 

197 


BBsO. 

00  197  lei. IS, i 
2FCYSCI) •GE*8B>  GOTO  299 
GOTO  299 
BBeYS(I) 

K*I 

CONTINJE 

CONTINJE  ^ 

PRINT*, “MAX  ON  OVER  CONDITION  #  "»BB 
PRINT*, “AT  SAMPLE  RATE  OF  “,TSSCK> 
CALL  PL0T(9*,-4.,-3) 

CALL  H3RAPH  (TSS,YS, IS, IYSS,  1,0,1) 

CALL  PLOTE(Ml) 

STOP 
END  . 


SUBROUTINE  C0N0(F2,G2,C,N,Z0IM,Q) 

DIMENSION  F2(I0IM,XDIP),G2(IDIM,IDIH) 

DIMENSION  C1C10,1C)  ,C 2<1 0,10 > ,C3(01, 10>  ,CM10,10> 
DIMENSION  C5(10, 10)  ,C  (IDIM.IDIM) 

DIMENSION  RC5(10«10) 

DIMENSION  P(10,10) 

DIMENSION  G7(10,10),Q(10),HK  (20)  ,B(10,  10) 

CALL  C0PY(C,CS,N,I0IM) 

CALL  COPY (G2,G7,N«X01M) 

C  FIND  HANKEL  C IMPULSE  RESPONSE)  MATRIX 
00  1115  J®1 ,N, 1 
C2(l, J) *C(1,J) 

CA(Jfl) *G7(J,1) 

1119  CONTINUE 

00  1113  1*2, N,  1 

CALL  MJLT (C5,F2, Cl, 1, N,N ,101 H) 

CALL  MJLT (F2«G7,C3,N,N,1,I0IM) 

00  111%  J*1,N, 1 
C2(I, J) *C1(1, J) 

CL(J,I)*C3(J,1) 

1114  CONTINUE 

CALL  C0PY(C1,CS,N,IUXM) 

CALL  C0eY (C3,G7, N,IDZN) 

1113  CONTINUE 

CALL  MULT  (C2,C4,RC9 ,N»N, N,IOIM) 

XAelDIM 
XOGT  *N 

CALL  COPY fRC9,P,N,X0XH) 

CALL  LSVDF(P,XOZM,NyN|B, •1,*1,Q,NX,IER) 

RETURN 

END 


101 


i 


C 

C 

SUBROUTINE  EFT  (A ,  B,  N,  102  M,A4,B5,M,  OELTI 
C  THIS  SUBROUTINE  FXNO  NY  OXSCRETE  SYSTEM  F  ANO  6 
C  FROM  A  ANO  3.  A  IS  NXN  ANO  B  XS  NXN 
C  DELT  IS  HY  TINE  INCREMENT •  ANO  M  IS  THE  NUMBER 
C  OF  ITERATIONS  IN  HY  SUM  I  WANT  TO  60* 

01 MENS I ON  B (IOXM , IOIM) ,A (I01H,I0IH) 

01 MENS I ON  A4(I0IM,X0IH) , 65 (XOIM»IDIM> 

01  MENS  I  ON  A2(10,1C)  ,A3(10,13> 

01  HENS  I  ON  AI(1Q,1Q) 

DIMENSION  A5(10,10)  ,AINV(10.10) 

OIHENSXON  84(10,10) 

DIMENSION  P(10,10)  > 

01 ME NS I ON  WKAREA (200) 

C  SET  UP  I0ENITY  MATRIX 

CALL  COPY  (A  ,P,N,  IOIM) 

lAsIOIM 

IOGT=N 

00  1002  1=1, IOIM, 1 
00  1003  J=1,N, 1 
AI(I,J)=0.0 
A4(I»J)*O*0 
1003  CONTINUE 

AI(I,I)=1.0 
1002  CONTINUE 

C  FINE  0  SUCH  THAT  F*Q*I  ANO  G=Q*AINV*B 
C6*l«9 
0ET1=1. 0 
00  1111  1=1, M,1 
DET1*DET1*0ELT 
C6=C6*I 
ABLE*0ET1/C6 

CALL  MULT (AI,A,A2,N,N,N, IOIM) 

CALL  COPY (A2,AI,N,IDIM) 

CALL  MULTXK(A2, ABLE, A3, N,N, IOIM) 

CALL  A03ING(A3,A4,A5, N,N,IDIM) 

CALL  COPY  (A5,A4,N, IOIM) 

1111  CONTINUE 
C  FINO  AINV 

CALL  LINV2F (A, N, IA, AINV, IOGT , WKAREA, IER) 
CALL  MULT (A4, AINV, B4,N,N,K, IOIM) 

CALL  MULT(B4,B,B5,N,N»1,ICIN> 

00  1001  1=1,  N,1  .  . 

A4(X,I)=A4(I,I)+1«0 
1001  CONTINUE 

CALL  COPY (A,P,N, IOIM) 

C  A4  IS  MY  F  ANO  B5  IS  MY  G 
.  RETURN 
ENO 
C 
C 
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SUBROUTINE  PREOX CT ( FI VC »  X»N« IOIH* Y) 

C  THIS  SUBROUTINE  OOES  THE  OUTPUT  PREOXCTXON 
C  PHASE  OF  OPOEC. 

DIMENSION  FKIOXNf IOIN) 

DIMENSION  xaOXM'XOXM)  ,Y  CIOIM,IOIN) 

DIMENSION  Zl<10#  10)  ,Z2<1 0»10) ,RC10 ,10) 

DIMENSION  C (IOIH* IOIH) 

1 

M7«N»N-1 

DO  196  II*N*N7fl 

CALL  MPONPCFlf II »  N» IOIMy R> 

CALL  HULT(RfX«ZlfN»N»  lflOIM) 

CALL  MULT(C»ZltZ2vlfNflfI0IM> 

Y(H4rl>-Z2(l»l) 

NL»N%«-1 

196  CONTINUE 
RETURN 
ENO 
C 
C 

C 

C 

SUBROUTINE  A00IN6 (A  »B >Ct Nf N« IOIH) 

C  N  IS  RON,  M  IS  COLUMN 
C  THIS  ADDS  TWO  MATRICES  OF  SAME  SIZE 

DIMENSION  A  (IOIH* IOIH)  »B  UOXN.IDIMJ  ,CtI3IM, IOIH) 

00  906  J*i»M» 1 
DO  906  I*1*N,1 
C(I»J)*MI|J)*B(X|J) 

906  CONTINUE 
RETURN 
ENO 
C 

c 

£«»+•» ♦*#******•***»***#•♦• **•»♦♦♦♦♦**** •******f *»*•**♦***•♦*** 

c 

c 

SUBROUTINE  COPY (A vB,N,XOXH> 

C  MUST  BE  A  SQUARE  MATRIX 
C  COPIES  A  INTO  B 

01 MENS ION  A (lOXMf XOIM) (B (XOXM»XQXNI 
00  1100  US*l»Nfl 
DO  1100  JT*1#N»1 
8(  JS,JT)  *A(  JS,  JT) 

1190  CONTINUE 
RETURN 
ENO 

6 

C 

6 

e 


SUBROUTINE  PRNMA (AyNy MylCIM) 

C  THIS  SUBROUTINE  PRINTS  OUT  ANT  SIZE  MATRIX 
C  H  IS  THE  NUMBER  OF  COLUMNS 
C  N  IS  THE  NUMBER  OF  ROMS 

OINENSION  A  (IOIHylOIN) 

00  1112  J*l»Nyl 
NRITE(Sylll) (A(Jyl) yl*ly Myl) 

111  FORMAT (*  Mf10(2XyE10.4)) 

1112  CONTINUE 
RETURN 
ENO 
C 
C 

C*******  ************************************  ***•****♦**♦•****•*<♦ 

c 

c 

SUBROUTINE  MPOHP (MyNPyNy lOIMyR) 

C  FINOS  R*M*  *NP 

01  MENS ION  M (IOIM y IOIH ) ,R (101 M, IDIM) 

01 MENS I ON  Rl(lOylO) yR2(l0yl0) 

00  193  J*lyIOIN,l 
00  19t  I*ly IOIMy 1 
R1(I,J)*0.0 
R2(I,J)*0.0 

194  CONTINUE 
RlCJyJMl.9 

193  CONTINUE 

00  195  JJ*l,NPyl 

CALL  MULT(H«Rl»R2yNyNyNyI0IM> 

CALL  COPT  (RZyRly  Ny  IOIM) 

195  CONTINUE 

CALL  COPT (R2*R,NyI0IM) 

RETURN 

ENO 

0 

c 

C 

c 

SUBROUTINE  NOIZE (RMSNOTSy OUTHEANyMN) 

C  SUBTOUTINE  NOISE  CALCULATES  THE  VALUES  OF  THE  MEASUREMENT  NOISE 
C  COMPONENTS  USING  A  RANDOM  NUMBER  GENERATOR  MODELLED  AS  SAUSSIAN 

GAUSS* 0* 

00  333  I*ly 12y 1 
CAUSS*GAUSS*RANF(Z10> 

933  CONTINUE 

GAUSS*  GAUSS  *5 • ♦OUTHEA N 
NN*GAUSS*RH SNOTS 
RETURN 
ENO 
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SUBROUTINE  HGRAPH(X,Y,N,  IC,NO»NP,NS) 

C  IF  IO(1)*O.OOC  BOX  IN  UPPER  RIGHT  CORNER 
C  IS  NOT  PLOTTED 

DIMENSION  X(l)  ,V(1)  ,XD(1)  f  IF(N0.E0.2)  GO  TO  30 
IF  (NO.LT.O)  GO  TO  1C 

CALL  SCALE(X,7.,N,1)  S  CALL  SC«LE( Y,C • ,N,1) 

10  CALL  PLOT  (8.5,0.  ,-3)  1  CALL  PLOT  (C  .  ,11.  ,3) 

CALL  PLOT  (-1.39,1.35,  3) 

CALL  PLOT (-7.15, 1.35, 2)  S  CALL  PL0T(-7.15,9.65,2I 

iF(iod) .en.cao)  go  to  25 

CALL  PLOT  (-7.05,  5.55,  3)  X  CALL  PLOT  (  -7. 05,7.55, 2) 

00  20  1*1, 7, 2 

20  CALL  SYMBOL  (I*. 1-6. 9,7  .85  ,.07,19(1)  ,90., 20 

CALL  PLOT  (-7.05,7.55,  3)  £  CALL  PL OT  (  «G.  05,7. 55,2) 
CALL  PLOT  (-6.09,9.55,2)  S  CALL  PLOT(-7.05,9.»S,2» 

CALL  PLOT (-7. 15,9. A5, 3) 

25  CALL  PLOT (-1.35, 9.65, 2)  S  CALL  PL0T(-1.35,1.35,2) 
CALL  SYM«OL  (-&.6r- ,1.1£  , .1,19  (13),  3. ,ti» 

CALL  AXIS  (-1.05, 2. 1,1 0(9)  ,-2J, 7., 93.  ,X  < N*1 )  ,  X(N  +  2)  ) 
CALL  AXIS(-l.fl5,  2.1, 10(11),  20,5.  .lSt  .,Y  (N+l)  ,Y(N*Z)I 
30  Y(N+2) *-Y (N+2) 

X(N+1)  *X  (N*l)  -2. 1*X  (N+2)  S  Y  (N+l)*  Y(N*-1  )*1.85*Y  (N+2) 
CALL  LINE(V,X,N, 1,NP,NS) 

X(N*1)  =X  (N+l)  *2,  1*X  (N*2)  t  Y(N41)*Y(N*1)-1.05*Y(N42I 
Y(N+2) *-Y (N  *2) 

RETURN  S  END 

SUBROUTINE  VGRAPH (X ,Y ,N, IC,NO,NP,NS) 

01 MENS I ON  X (1) ,Y (1) ,10(1 )  S  IP(NO. E9.2)  GC  TO  30 
IF  (NO.LT.O)  GO  TO  10 

CALL  SCALF(X,4.9,N,1)  1  CALL  SCALE  (Y,7 •  G»N,  1) 

10  CALL  PLOT  (8  .5, 0.  ,  -3)  i  CALL  PLOT  (0. ,11. ,3) 

CALL  PLOT (-1.35, 1.35, 3) 

CALL  PLOT  (-7. 15, 1.35,  2)  f  CALL  PLOT(-%  15 ,9.65 ,2) 
CALL  PLOT  (-1.35,9.155,  2)  1  IF  (10(1)  .  EO.  0  OT)  GOTO  25 
CALL  PLOT (-1.45,9.95, 3)  1  CALL  PL0T(-3. 45,9.55,2) 

00  20  I*  1,7,2 

20  CALL  SYMBOL  (-3.1  f  ,5  .4-1*  •  1G*  .07,10  (I) ,  9  .,20) 

CALL  PLOT  (-3.45, 9. 55,  3)  i  CALL  FL0T(-3.45, 8.55,2) 

CALL  PLOT (-1.45, 8.55,2)  *  CALL  PL0T(-1.A5,9.55,2) 

CALL  PLOT (-1.35, 9.6b, 3) 

25  CALL  PLOT (-1.35, 1 *3f , 2) 

CALL  SYMBOL (-6.65, 1.15, .1,10 (13), 0., A3) 

CALL  AXIS (-6.**,  1.85, 1 0(9  )  ,-20,4.3,  I.  ,X(N+1)  ,X(N*2)) 
CALL  AXIS(-G.4, 1.85,10(11),  20, 7.0*  90. ,Y  (N*l)  ,Y(N*2)  I 
SO  X(N4l)«X(N41)«-6.4*X(N42>  $  Y(N+1)  «Y  (MM) -1.85»Y (N*2> 

CALL  LTNE(X,Y9N,1VNP,NS) 

X(N+1) *X (N+l) -b. ***X  (N  +  2)  t  YCN+1)  «Y(N+1)41.85*V(NR2) 
RETURN  1  EMO 
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Appendix  B_ 


OPPEC  Program 

This  Appendix  contains  the  computer  listing  and  a  users  guide  for 
the  impllmentation  of  OPDEC.  The  program  reads  in  user  supplied  inputs 
and  generates  the  initial  states,  x(0).  The  next  phase  is  the  output 
prediction  phase.  In  this  phase  the  current  system  state,  x(k),  is  used 
to  predict  the  outputs  in  the  future.  This  predicted  output  is  then 
used  to  calculate  the  next  input,  using  OPDEC 's  control  lav.  The  input 
is  then  applied  to  the  true  system  to  update  the  state  vector  one 
sample  step.  This  updated  state  vector  is  then  used  in  the  prediction 
phase  on  the  next  cycle  through  the  program. 

Program  Inputs  and  Initial  Conditions 

The  user  must  supply  the  program  with  the  following  data.  All  data, 
unless  specified  otherwise,  is  read  in  using  an  unformated  read  statement. 

1.  A  title  to  be  used  on  the  output  plot.  This  title  is  usually 
the  sample  rate  the  system  is  using.  It  is  used  to  identify  output 
plots.  The  title  is  limited  to  twenty  spaces  and  read  in  using  an 
alphanumeric  format. 

2.  The  program  output  contains  two  plots.  One  plot  is  of  the 
system  output  and  the  other  is  of  the  inputs  calculated.  Each 
plot  has  a  title  box  in  the  upper  right  hand  corner.  This  program 
has  the  option  of  drawing  this  title  box  or  not.  The  next  input, 
therefore,  is  a  switch  variable.  If  this  variable  is  less  than  or 
equal  to  sero  both  title  boxes  will  be  drawn.  If  it  is  less  than 
one  hut  greater  than  tero  the  title  box  will  be  drewn  only  on  the 
control  input  plot  and  not  on  the  output  plot.  If  the  input  is 
greeter  than  one,  no  title  box  will  be  drawn  on  either  plot. 


3.  System  size  n.  n  is  an  integer  value  and  must  be  less  than  or 
equal  to  ten. 

4.  A,  B,  and  C  from  the  true  system  state  matrix  equation 

x  -  A  x  +  B  u 
y(t)  -  C  x(t) 

A,  B,  C  are  all  read  in  as  n  by  n  matrices  and  must  be  discretizable 
at  the  users  supplied  sample  rate.  A  must  also  be  invertible. 

5.  Am,  Bm,  and  Cm  from  the  perturbed  system  state  matrix  equation 

x  ■  Am  x  +  Bo  u 
y(t)  -  Cm  x(t) 

Am,  Bm,  Cm  are  all  read  in  as  n  by  n  matrices  and  must  be  discret¬ 
izable  at  the  users  supplied  sample  rate.  Am  must  be  invertible. 

6.  Sample  rate  DELT,  which  the  system  is  to  use.  DELT  is  a  real 
number. 

7.  The  program  has  an  option  of  using  a  smaller  model  than  the 
system  order,  n,  in  calculation  of  the  next  input.  This  option 
is  the  next  input  to  the  program,  N6.  This  tells  the  program  the 
size  of  the  Hankel  matrix  to  use  in  the  calculation  of  the  closed 
loop  input.  M6  is  an  Intiger  and  must  be  less  than  or  equal  to 
the  system  size,  n,  and  greater  than  or  equal  to  1. 

8.  The  next  input  is  an  Intiger,  and  tells  the  program  the  number 
of  cycles  to  implement  OPDEC's  control  law. 

9.  The  next  input  is  the  standard  deviation  of  the  noise  one  wants 
to  add  to  all  the  states. 

10.  US  la  the  standard  deviations  of  the  guaaslan  noise  added  to 
the  input .  If  this  value  is  leas  than  or  equal  to  aero,  then  no 
input  noise  Is  added. 
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11.  The  next  inputs  are  three  logic  switches  used  in  making  the 
program  follow  a  path  the  user  desires.  The  inputs  are  read  in 
using  a  logical  format.  The  logic  switches  can  take  on  two  values: 
a  T  (true)  or  an  F  (false).  The  first  logic  switch,  SW1,  if  true 
tells  the  program  to  use  the  true  model  in  the  prediction  phase 
of  the  program.  If  SW1  is  false,  the  program  uses  the  perturbed 
model  in  the  prediction  phase  of  the  program.  The  second  logic 
switch,  SW2,  if  true  tells  the  program  to  use  the  true  model  in 
the  control  phase  of  the  program.  If  SW2  is  false  the  program 
uses  the  perturbed  model  in  the  control  phase  of  the  program. 

SW3  is  the  third  logic  switch.  If  true  it  will  add  random  state 
noise  with  a  standard  deviation  selected  by  the  user  (Input  9) . 

If  SW3  is  false  no  state  noise  will  be  added. 

After  reading  in  the  inputs,  the  program  discretizes  the  true 
system  and  the  perturbed  system  using  subroutine  EFT.  The  next  sections 
are  done  sequenclally  in  the  closed  loop  control  loop. 

Prediction  Phase 

This  phase  uses  SW1  to  select  which  model  to  use  in  the  output 
prediction.  The  output  prediction  is  done  by  using  subroutine  PREDICT. 
Also  done  in  this  phase  is  to  set  up  some  titles  to  be  drawn.  These 
title  changes  inform  the  user  which  model  was  used  in  the  prediction 
phase. 

Control  Phase 

This  phase  uses  SW2  to  select  which  model  to  use  in  the  control 
phase.  The  output  prediction  is  used  to  determine  an  input.  This  is 
done  by  using  subroutine  CONTR.  There  are  also  some  additional  title 
changes,  for  plotting  purposes,  to  Inform  the  user  on  the  output  plots 


108 


■s*' 


m 


avia sewiMHier  ^ 


which  model  was  used  in  the  control  phase. 


System  Implementation 

Uhat  is  done  in  this  phase  is  to  take  the  input  calculated,  state 
vector  and  apply  them  to  the  true  system  to  update  the  state  vector  one 
sample  step.  This  is  done  by  using  subroutine  TRUTH.  Also  in  this 
phase  SW3  is  used  to  add  state  noise  if  the  user  wants  state  noise  added. 
There  are  some  additional  title  changes  to  inform  the  user  on  the  output 
plots  if  state  noise,  input  noise  or  both  was  added. 

Output  Plots 

After  the  closed  loop  control  is  done  the  system  then  plots  the 
sampled  output  versus  time  and  control  inputs  versus  time.  The  plotting 
is  done  using  subroutine  HGRAPH  or  VGRAPH. 

Major  Subroutines 

Subroutine  CONTR(F2.G2.C.Y.N.IDIM,D,N6 ,KZ ,RC15) 

This  subroutine  uses  the  supplied  discrete  system  F2,  G2,  C  and 
creates  the  Hankel  matrix.  The  Hankel  matrix  is  then  Inverted  and  is 
output  in  RC15.  This  matrix  is  then  multiplied  by  the  output  prediction 
vector,  Y,  to  determine  the  input  D.  Subroutine  CONTR  is  called  many 
times  during  the  closed  loop  control  phase.  So  to  save  computer  time 
the  Hankel  matrix  is  created  and  inverted  the  first  time  CONTR  is  called 
and  then  stored  in  memory,  so  in  subsequent  callings  the  matrix  already 
exists  and  does  not  need  to  be  recomputed.  N6  was  described  in  the 
input  section  of  this  appendix.  IDIM  is  the  initial  dimension  of 
F2,  G2,  C  and  RC15. 

Subroutine  PREDICTfFl.C.X.H.IDIM.Y) 

This  subroutine  uses  parts  of  the  discrete  system  model,  FI,  C 
and  state  vector  x  to  predict  the  output  at  discrete  time  n  to  2n-l. 


( 


( 


( 


These  outputs  are  then  put  into  vector  Y.  IDIM  is  the  initial  dimen¬ 
sion  of  FI,  C,  X,  and  Y. 


Subroutine  EFT(A,B,N, IDIM, F,G,M, PELT) 

This  subroutine  finds  the  discrete  system  F  and  G,  from  a,  n  by  n, 
and  B,  n  by  n,  at  the  sample  rate  of  DELT. 


T  4-  ?  a*  (DELT)j 

1  i-i  Tiyr 


where  M  is  the  number  of  sumations  the  user  wants.  Then  G  is  calculated 
via 

G  **  (F  -  I)  A'1  B 

where  I  is  the  identity  matrix.  This  equation  is  the  reason  why  the 
restriction  that  A  and  Am  must  be  invertable  was  stated  in  the  input 
section  of  this  appendix. 

Subroutine  MPOWP (M, NP , N , IDIM, R) 

This  subroutine  takes  matrix  M,  n  by  n,  to  the  power  of  NP.  Np 
is  an  integer  value.  The  answer  is  then  put  into  matrix  R.  IDIM  is 
the  Initial  dimension  of  M  and  R. 

Subroutine  TRUTH(A.B.C.DELT,D,X,IDIM,N,TSS,YS,IS,US,RMS) 

This  subroutine  takes  the  input  D  and  state  vector  x  and  applies 
them  to  the  true  system,  A,  B,  C  to  update  the  state  vector  one  sample 
rate.  To  insure  the  complete  output  response  and  help  smooth  output 
plots  the  true  system  runs  at  a  sample  rate  that  is  eleven  times  faster 
than  the  user's  selected  sample  rate.  First  the  true  system  is 
discretized  using  subroutine  EFT,  at  a  sample  rate  of  DELL,  where 

DELL  -  DELT/11.0 

Next,  input  noise  is  added  if  RMS  is  greater  than  zero.  Then  a  white 
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guassian  noise  will  be  added  to  the  input.  The  noise  will  have  a 

standaru  deviation  of  RMS.  Next  the  system  states  are  updated  one 

sample  rate.  Also  the  arrays  to  be  plotted  TSS,  US,  and  YS  are 

created.  IS  is  a  counter  for  setting  up  these  arrays.  Then  after 
* 

the  state  vector  has  been  updated  the  time,  state  vector  and  output 
are  printed. 

Subroutine  NOIZE (RMSNOTS , OUTMEAN , WN) 

This  subroutine  uses  a  random  number  generator  to  create  a  white 
guassian  noise  with  a  standard  deviation  of  RMSNOTS  and  a  mean  of 
OUTMEAN.  The  noise  is  output  in  WN. 

Closing  Remarks 

This  program  can  only  be  used  with  a  single  input,  single  output 
system.  But  all  of  the  subroutines  and  the  vital  parts  in  the  main 
program  are  set  up  for  a  system  with  n  inputs  and  n  outputs.  The 
plotting  routines  are  set  up  so  that  with  minor  modifications  one  can 
have  multi-plots  per  run.  The  program  uses  some  subroutines  from 
IMSL  library  packages,  they  are  LINV2F  and  LSVDF. 


Ill 


PROGRAM  THESIS (INPUT, OUTPUT, TAPE5=INPUT,TAPE6*OUTPUr, PLOT) 
DIMENSION  AT(10,10) ,3T(1C,1G> ,CT(i0,10) 

DIMENSION  AM(13 , 1C) ,BM (1 C  ,10 ) , CM  CIO, 10 ) 

DIMENSION  X0C10, 10), Y(10, 10) 

DIMENSION  X(13,10),RC15(1G,10) 

C  MAX  NUMBER  OF  INPUTS  ONE  CAN  CALCULATE  WITH  DIMENSION  Or  1060 
C  IS  58  BECAUSE  OF  OIVISION  BY  11  IN  TRUTH  SUBROUTINE 
DIMENSION  TSS(lCCO)  ,YS(1000)  , US (10 00) 

01  HENS I ON  IYSS(17),I0SS(17)«I?SS(17),IXSS(17),ISSS(17) 
DIMENSION  CON 1(10, 10) ,CCN2( 10,10) 

DIMENSION  FT(10, 10) ,FM(1 0 ,10 ) ,GT (1 0, 10) ,GK(10,10) 

COMHON/MA IN1/NDIH ,NCIM1, C0N1/MAIN2/C0M2/IN0U/KIN,K0JT ,< PUNCH 
LOGICAL  SW1,SW2,SW3 
C  SET  UP  TITLES  FOR  CALCOMP  PLOTS 

DATA  IYSS  (D/2QH  SAMPLE 0  OUTPUT  / 

OATA  IYSS (3)/20H  TIME  IN  SECONDS  / 

OATA  I  YSS  ( 1 1)  /  23  H  SAMPLED  OUTPUT  / 

OATA  IOSS ( 1 ) /20H  CONTROL  INPUTS  / 

OATA  IOSS (9) /20H  TIME  IN  SECONDS  / 

OATA  IOSS (11) /2GH  CONTROL  INPUTS  / 

OATA  IZSS (D/20H  TRUE  PREDICTION  / 

OATA  I7SS(3)/20H  FALSE  PREDICTION  / 

OATA  IZSS (5)/20HTRUE  CONTROL  MODEL  / 

OATA  I73S(7)/2QHFALSE  CONTROL  MODEL  / 

DATA  IZSS(9)/40HSAMPLED  OUTPUT  WITH  STATE  NOISE  ADDED  / 
DATA  IZSS (13) /49HSA HP LEO  OUTPUT  WITH  NO  NOISE  ADDED  / 

OATA  IXSS (1) /40HCONTROL  INPUTS  WITH  STATE  NOISE  ADDED  / 

OATA  IXSS (5)/L JHCONTROL  INPUTS  WITH  NO  NOISE  ADOED  / 

DATA  IXSS (9)/40HOUTPUT  WITH  INPUT  AND  STATE  NOISE  A3 DEO  / 
OATA  IXSS  (13)  /40HCONTROL  INPUT  WITH  IN°UT  AND  STATE  NOISE/ 
OATA  ISSS  (1)/40H SAMPLED  OUTPUT  WITH  IN=»UT  NOISE  ADDED  / 
OATA  ISSS (5 )/% OH CONTROL  INPUT  WITH  INPUT  NOISE  ADDED  / 

C  READING  IN  4TH  LINE  OF  UPPER  RIGHT  HAND  BOX  OF  PLOTS 
C  USUALLY  READ  IN  IS  "SAMPLE  KATE  OF  XXXXX" 

RE AD (5, 5 5£) IYSS ( 3) , IYSS( 4) 

556  FORMAT (IX, 2A10) 

IOSS (3) SI YSS (3) 

IOSS (4 ) si YSS (4) 

C  RE  AO  A  NUMBER  TO  HAKE  BOX  IN  UPPER  LEFT  HAMO  CORNER 
C  OR  NOT  I  ZZZ.GT.1.0  THEN  NO  BOX  IN  CONTROL 
C  IF  ZZZ.GT.3,  BUT.LE.i.O  THEN  NO  BOX  IN  OUTPUT  PLOT* 

C  BUT  SOX  IN  CONTROL  PLOT 
REA0*,7Z7 

IF(ZZZ.GT.0.0)IYSS(1)=0.0 
IF(ZZZ*GT*1.0)IDSS(1)=0*0 
C  READING  IN  SIZE  OF  MV  SYSTEM  N  ' 

C  N  MUST  BE  LESS  THAN  OR  EQUAL  TO  10 
READ*, N 
IDIM»10 
19*55 
NOIM*10 
NOIM1-11 


KIN=5 
KOUT  *5 
KPUNCH=Z 

C  REAOZNG  IN  TRUE  SYSTEM  IN  STATE  MATRIX  FORM 
REAO*,  <(AT(I,J) , J=1,N) ,1  =  1,N ) 

READ*, <(BT(I,J) , J=1,N) ,I=i,N> 

REAO*, ((OT(I,J),J=l,N),I=l,N) 

C  REAO  IN  PERTURBED  SYSTEM  IN  STATE  MATRIX  FORM 
READ'S  ((AM(I,J),J=1,N),I=1,N) 

READ*, ((RM(IfJ),J=i,N)#Isi*N) 

READ** <(CM(I,J), J=1,N) ,1=1, N> 

12  CONTINUE 
C  REAO  IN  SAMPLE  RATE  USEO 
READ*, DELT 

C  TO  STOP  PROGRAM  LET  DELT  BE  *LE.0*D 
IF(DELT.LE.O.O)  GOTO  96 

C  READING  SIZE  OF  HANKEL  MATRIX  USEO  IN  INVERSE 
C  IN  SUBROUTINE  CONTR 
REAO*, NS 

C  READING  IN  THE  NUMBER  OF  TIMES  YOU  CALCULATE  AN  INPUT 
REAO*,KZ 

PRINT*, “OELT=“,OELT 
C  REAO ING  RMS  VALUES  FOR  NOISES 
C  FIRST  IS  FOR  NOISE  AOOEO  TO  STATES 
C  THEN  NOISE  APOED  TO  INPUTS 

C  IF  RMS  LESS  THAN  ZERO  NO  NOISE  AOOEO  TO  INPUTS 
READ*, RMS NO  IS, RMS 

C  REAO  IN  LOGIC  SMITH  FOR  COMPUTATION 

C  SW1  FOR  PREDICTION  SW1=TRUE  USING  MATRIX  AT  FOR  PREDICTION 
C  SW1= FALSE,  USING  MATRIX  AM  FOR  PREDICTION 
C  SW2= TRUE,  USING  TRUE  MOOEL  FOR  CONTROL  PART  OF  PROBLEM 
C  SN2*  FALSE,  USING  PERTURBEO  MODEL  FOR  CONTROL  PART  OF  PROBLEM 
C  $H3= TRUE,  USING  NOISE  IN  STATE  UPDATE 
C  SW3* FALSE,  NO  NOISE  BEING  USEO* 

REAO (5 , 555)  SM1,SH2,SH3 
555  FORMAT  (  3L1) 

C  SET  COUNTER  TO  ZERO 
IS=0 

C  GENERATE  MY  DISCRETE  F  ANO  G  MATRIX  FOR  TRUTH  AND  MODEL  SYSTEM 
Ml =5  5 

CALL  EFT  CAT ,BT,N, IOIM ,FT ,GT«Ml,OELT) 

CALL  EFTCAM,BM,N, IOIM, FM,GM, Ml, DELT) 

C  GENERATE  INITIAL  CONDITIONS  RANDOMLY 
CALL  RANSET <Z9> 

00  444  1=1, N, 1 
XO (I , 1) = 1C*RANF( Z9> 

444  CONTINUE 

PRINT*,*  * 

PRINT*,"  INITIAL  CONDITIONS  ARE" 

CALL  PRNMA (XO,N, 1,I0IM) 

CALL  COPY <XO,X,N, IOIM) 

PRINT*,"  " 


00  445  KNI*1,KZ,1 

C  GOTO  PROPER  PREDICTION  USING  LOGIC  SWITCH  #1 
IF (SWi) SOTO  29 
IYSS (5) *IZSS(3) 

IYSS(5>*IZSS<4> 

I0SS(5)=I7SS(3) 

I0SS(6)=I7SS(4) 

PRINT*, "USING  MOOEL  MATFIX  SW1  IS  FALSE" 

CALL  PREDICT (FM, CH,X,N,IOIH,Y) 

GOTO  30 

29  PRINT*, "USING  TRUTH  MATRIX  SWi  IS  TRUE" 

IVSS(5)=I7SS(1) 

lYSS(S) *I7SS(2) 

I0SS(5)=IZSS(i) 

I0SS(6)=I7SS(2) 

CALL  PREDICT (FT , CT, X, N,I OIM, Y) 

30  PRINT*,"  MY  PREDICTION  IS  " 

-CALL  P.RNMA  ( Y ,  N ,  1 , 1  DIM ) 

PRINT*,"  " 

C  GOTO  PROPER  CONTROL  USING  LOGIC  SWITCH  *2 
IFCSW2)  GOTO  39 
IYSS (7  ) -IZSS(7) 

IY SS(B)=IZSS(3) 

IDSS(7)=IZSS(7) 

ZOSS (8) -I^SS (8) 

PRINT*, "USING  HOOEL  MATRIX  FOR  CONTROL  SW2  IS  FALSE" 
CALL  CONTR(FM,GM,CM,Y ,N, IOIM,D,NS, KNI,RC15) 

GOTO  40 

39  PRINT*, "USING  TRUTH  MATRIX  FOR  CONTROL  SW2  IS  TRUE" 
IYSS (7 ) =1 ZSS(5) 

IYSS(8)=I7SS(6) 

IDSS (7) *I7SS(5) 

IDSS(8) =IZSS(6) 

CALL  C0NTR(FT ,GT ,CT ,Y ,N, I0IM,D,N6, KNI ,RC1F) 

40  PRINT*,"  " 

PRINT*, "MY  INPUT  IS  ",0 
PRINT*,"  " 

CALL  TR'JTH(AT,3T  ,CT  ,0ELT  ,O,X0,IOIM,N,TSS,YS,IS,US,R1S) 
C  PUT  IN  NOISE  USING  SWITCH  #3 
IF(RMS.LE.O.)  GOTO  798 
IF (SW3)  GOTO  49 
PRINT*,"  " 

PRINT*, "INPUT  NOISE  ADOEO  BUT  NO  STATE  NOISE  AOOED" 

00  86S  KN=1 3,16, 1 
IYSS (KN) *ISSS (KN-12) 

I0SS (KN) *ISSS (KN*8) 

066  CONTINUE 

CALL  COPY(XO»X,N,IDIN) 

GOTO  445 

49  PRINT*,"  " 

PRINT*, "BOTH  INPUT  NOISE  ANO  STATE  NOISE  AOOED" 

00  067  KN»13,16,1 
IYSS (KN) >IXSS (KN*4) 
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IDSS (KM) *IXSS(KN) 

$67  CONTINUE 
GOTO  m 

79$  IF  (SMS)  GOTO  799 

PRINT*, "NO  INPUT  NOISE  OR  STATE  NOISE  ADDED" 

00  6 78  <N=13,16, 1 

IYSS«N)*IZSS<KN) 

IOSS(KN)  *IXSS(KN-8) 

$7$  CONTINUE 

CALL  COPY (XO,X,N,IOIH) 

GOTO  44$ 

799  PRINT*, "NO  INPUT  NOISE  BUT  STATE  NOISE  AODEO" 

00  $7$  KN=13,16, 1 
IYSS(KN)*IZSS<KN-4) 

IOSS (KN) *IXSS( KN*12) 

$79  CONTINUE 
77$  CONTINUE 

DO  446  KP=i,N,l 
OUTHEANsO.O 

CALL  NOIZE(RHSNOIS,OUTNEAN,HN) 
X(KP,1)=XC(KP,1)  ♦HN 
446  CONTINUE 

PRINT*, "STATES  AFTER  NOISE  ADDED" 

CALL  PRNMA (X,N,1, IOIH) 

445  PRINT*,"  " 

CALL  PL0T($.,-4.,~3> 

CALL  HGR« PH(TSS,YS, IS, IYSS, 1,0,1) 

CALL  PLOT <0«, -4,, -3) 

CALL  HGRAPH(TSS,US,IS,IDSS,1,0,1) 

GOTO  12 

96  CONTINUE 

CALL  PL0TECM1) 

STOP 

END 
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SUBROUTINE  CONTR <F2,G2,C , Y,N .IOIM,  0,  N5 ,K?,RC15) 

C  THIS  SUBROUTINE  CREATES  THE  HANKEL  HATRIX 
C  ANO  USES  IT  FOR  COMPUTING  THE  INPUT  MEEOEO 
DIMENSION  F2(I0IM,I0IM),G2(  IOIM,  IOIM) 

DIMENSION  CldOi  1C)  ,C2(10,10)  ,C3(01,10)  ,C*  (10,10) 
DIMENSION  C5(10,10)  ,C  (IOIP,IDIM) 

DIMENSION  RC5( 10 , 10) , RC1 5 (IDIM,IOIM) 

DIMENSION  P(lOflO) 

DIMENSION  G7(10,10),0(1G) ,WK (20) ,8(10,10) 

DIMENSION  Y (IOIM , IOIM) ,U  (10, 10) 

IF(KZ.GT.l)  GOTO  666 
CALL  COPY (C,C5,N,I0IM) 

CALL  COPY (G2,G7,N,I0I K) 

C  FINO  HANKEL  (IMPULSE  RESPONSE)  MATRIX 
DO  1115  J=1,N, 1 
C2 (1, J) =C (1 , J) 

C<t  (J,l)-G7  ( J,l) 

1115  CONTINUE 

00  1113  1-2, N,1 

CALL  MULT (C5*F2, Cl, 1, N,N »IDIM) 

CALL  MULT (F2,G7,C3,N,N,1 , IOIM) 

DO  1114  J=1 ,N, 1 
C2  (I  ,  J)  =  C1  ( 1,  J) 

C4(J,I)  *C3(  J,l) 

1114  CONTINUE 

CALL  COPY  (C1,C5,  N,IDIH) 

CALL  COPY (C3,G7 , N,IOIM) 

1113  CONTINUE 

CALL  MULT (C2,C4,RC5,N,N, N,IDIM) 

IA*IOIH 

IOGT*N 

PRINT*, "HANKEL  MATRIX  IS" 

CALL  PRNMA (RC5,N6,N6, IOIM) 

CALL  COPY (RC5,P,N, IOIM) 

CALL  LSVDF(P,IOIH,N,N,B,-ly-l,Q,HK,IER) 

PRINT*,"  " 

PRINT*,"  SINGULAR  VALUES  OF  HANKEL  MATRIX" 

CALL  PRNMA(Q,N,l,IOIM) 

CALL  GMINV(N6,N6,RC5,RC15,MR,5) 

PRINT*, "RANK  OF  HANKEL  MATRIX  IS  ",  MR 
*6*  CONTINUE 

CALL  MULT (RC15,Y,U,N6,N6 ,1,10IM) 

0>-U(N6,l) 

RETURN 

END 

C 

e 

£+***»*»*»»»»»»»****«**+#**«*»****•**»»»*•*»»«*»••««»•«*«»»* 
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SUBROUTINE  PREOICT(Fl,Cy XyN«IDXHyY> 

C  THIS  SUBROUTINE  DOES  THE  OUTPUT  PREDICTION 
C  PHASE  OF  OPOEC. 

OX HENSION  Fl(IOIHylOIK) 

01 MENS ION  X(IOIH,IOIH) ,Y  <IOIH,IOIN> 

01 HENS I ON  Zl(10f 10) yZZUGylO)  yR(lOylO) 

01 HENS I ON  C (IOIH, IOIH) 

Nbtl 

N7*N»N-i 

00  196  IX*NyN7,l 

CALL  HPO«P(FltH,N,  IOIMyR) 

CALL  MULT(RyX,Zi,N,N, 1,I0IH) 

CALL  HULT (CyZlyZ2,lyN,ly IOIH) 

T(N4rl)=Z2(l,l) 

N6*N4+1 
196  CONTINUE 
RETURN 
END 
C 
C 

C 

C 

SUBROUTINE  AOOING (A  ,B  ,C  y  N  ,H« IOIH)  j 
C  N  IS  ROH,  N  IS  COLUMN 
C  THIS  ADOS  TWO  MATRICES  OF  SAME  SIZE 

01 HENSION  A  (IOIH, IOIH  >,B(IDIM,IOIM),C(IOIK,  IOIH) 
DO  906  J*1 y  H, 1 
DO  906  1*1, N,1 
C(I,J)*MI,J)»B(I,J) 

906  CONTINUE 
RETURN 
END 
C 
C 

£»**«» *******,**»*»**•*#♦*# •**»*«***»»**»«»** 4* *«*«*«« 

c 

c 

SUBROUTINE  COPY( AyByNyIDIH) 

C  MUST  9E  A  SQUARE  MATRIX 
C  COPIES  A  INTO  8 

DIMENSION  A (IOIH , IOIH ) ,B  (IOIH, IOIH) 

00  1100  JS*lyNrl 
00  1100  JT*lyNl,y  1 
8<  JSyJT)  *A(  JSfUT) 

1188  CONTINUE  "  / 

RETURN 

END 


c 

c 

SUBROUTINE  EFT(A,B,N,X01K,A4,B5,M,0ELT) 

C  THIS  SUBROUTINE  FIND  HY  DISCRETE  SYSTEM  F  AND  6 
C  FROM  A  ANO  9.  A  IS  NXN  ANO  B  IS  NXN 
C  CELT  IS  HY  TIHE  INCREMENT,  ANO  H  IS  THE  NUHBER 
C  OF  ITERATIONS  IN  HY  SUH  I  WANT  TO  60* 

01  HENS  I  ON  B(XOXM,IOIM)  ,A  (IOIN,IDIH> 

DIMENSION  A4(I0IM,IDI H) , 65 (I DIM, X DIM) 

01  HENS  I  ON  A2(10,1C)  ,A3(10,13) 

OXHENSION  AK1Q,  10) 

OXHENSION  A5(10, 10) ,AINV (10, 10) 

01 HENS ION  84(10,10) 

OXHENSION  P(10,10)  > 

OXHENSION  MKAREA (200) 

C  SET  UP  IOENITY  MATRIX 

CALL  COPY(A  ,P,N,  XOXH) 

XA*IDIM 

IDGT*N 

00  1002  I*i,IOIM,l 
00  1003  J*1,N, 1 
AI(I,J)=0,0 
A4(I,J)*0.0 
1003  CONTINUE 

AI(IfI)-1.0 
1002  CONTINUE 

C  FINE  Q  SUCH  THAT  F*Q*X  ANO  G*Q*AINV*8 
C6*l«3 
OET 1=1.0 
00  1111  1*1 «M, 1 
DET1*DET1*0ELT 
C6=C6*I 
ABLE*DET1/C6 

CALL  HULT (AI,A,A2,N,N,N, XOXH) 

CALL  COPY (A2,AI,N,I0IH) 

CALL  MULTXK(A2,ABLE,A3,N,N,IDIM) 

CALL  AD3ING(A3,A4,A5,N,N,IDIM) 

CALL  COPY (AS, A4, N, XOXH) 

1111  CONTINUE 
C  FIND  AINV 

CALL  LXNV2F (A,N, I A, AINV, IDGT ,WKAREA, IER) 

CALL  HULT (A4,AXNV,B4, N,N,N,IOXH) 

CALL  HULT (B4,B,B5,N,N ,1, XCXH) 

00  1001  1*1,  N,1 
A4(I,I)*A4(I,I)+1«0 
1001  CONTINUE 

CALL  COPY  (A,P,N,ZOIN) 

C  A4  IS  MY  F  ANO  B5  IS  MY  G 
.  RETURN 
END 
C 

c 

c*********** •••****•*••***♦***♦♦***♦**•********• •*••♦#♦****•♦• 
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c 

SUBROUTINE  TRANSP (A ,B , N, IDIM) 

C  THIS  TRANSPOSES  A  ANO  PUTS  INTO  B 
OIHENSION  A  (IDIM, IOIH) 

DIMENSION  B  (IDIM, IDIM) 

00  13CC  J=1  ,N, 1 
DO  1309  1*1,  N,1 
Bdy  J)  *A(  J,I) 

1300  CONTINUE 
RETURN 
END 
C 
C 

C 

C 

SUBROUTINE  MULTXK(A,D,C,M,N,IOIM) 

C  THIS  MULTPLIFS  A  MATRIX  BY  A  CONSTANT 

DIMENSION  A UDIM, IOIH), C (IOIH, IDIM)  . 

DO  90S  J*1,M,1 
00  905  1=1, N,1 
C( J,  I)  *A  (  J, I) *0 
905  CONTINUE 
RETURN 
END 
C 
C 

c 

c 

SUBROUTINE  MULT CS, ST, H,l ,H,N, IDIM) 

C  THIS  MULTPLIES  TWO  MATRICIES  TOGETHER 

OIHENSION  SdOIM, IDIM), ST(IOIM, IDIM)  ,H(  IDIM, IOIH) 

00  2000  1=1, 1,1 
00  2000  K=1,N,1 
SUM* 0,0 

00  2001  J=1,M,1 
SUM*SUM+S  d,  J)  *ST  (J,K) 

2001  CONTINUE 

H(I,K) *SUH 
2000  CONTINUE 
RETURN 
ENO 
C 

c 

£**«*»***»*»»»»*»*****»***»*•* *«*♦**»»»****»* «*»«*«»»**»*»*»» 

c 


c 

c 

c 

C  ,  "  '  /  ^  Z'  ,  ✓  /  ' 

SUBROUTINE  TRUTH (  Ai,Bl ,C8 , CELT, 0,X , IOIM, N, TSS, VS, IS,  US,  RMS) 
C  THIS  IS  THE  TRUE  SYSTEM  RESPONSE  WHICH  JS  RUNNING 
C  ELEVEN  TIMES  FASTR,  SO  THAT  THE  OUTPUT  PLOTS  ARE  SMOOTH 
C  IT  ALSO  CREATES  THE  ARRAYS  FOR  CALCOMP  PLOTS 
DIMENSION  X  (IDIM, I DIM) ,TSS(1C00) 

DIMENSION  A  1(10, 10) ,91(10,10) ,C8 (10,10) 

DIMENSION  26(10,10) ,27(10,10) 

DIMENSION  F3(10«10),G3(10,10) 

DIMENSION  SY(10, 10) ,YS(10C0) ,US(1000) 

M*50 

OELL*DELT/11.0 

CALL  EFT (A1,61,N,I0IM,F3 ,G3«H,DELL) 

00  702  IP=1,11,1 
IS=IS*1 

IF  (IS*  EQ*  1)  GOTO  6% 

TSS(IS) *TSS (IS-1) ♦OELL 
GOTO  kS 

hk  TSS  (IS)  =OELL 

45  CONTINUE 

e=D  .  •  •  - . 

IF(RMS.LE.O.O)  GOTO  777 
OUTM£AN=0* 

CALL  N0I7E(RMS,0UTMEAN,WN) 

BaD+HN 

777  CONTINUE 
US (IS) =8 
776  CONTINUE 

CALL  HULT(Fi3,X,26,N,N,l,I0IN> 

CALL  MULTXK  (G3 ,B, 27 ,N ,1, IDIM) 

CALL  A00ING(Z6,Z7,X,N,1, IOIM) 

CALL  MULT (C6,X,SY,1,N,1, IOIM) 

YS(IS) «SY (1,1) 

702  CONTINUE 

PRINT*,"  TIME  IS  **,T$S(IS) ,"  NY  STATES  ARE  " 

CALL  PRNNA(X,N,l,IOIM) 

PRINT*,"  " 

PRINT* ,"  OUTPUT  IS  ",YS(IS) 

PRINT*, "  " 

RETURN 

END 

C 

C 
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SUBROUTINE  PRNMA (A, N, M,1CIM) 

C  THIS  SUBROUTINE  PRINTS  OUT  ANY  SIZE  MATRIX 
C  H  IS  THE  NUMBER  OF  COLUMNS 
C  N  IS  THE  NUMBER  OF  ROMS 

OINENSIDN  A  (IOIH« IOIM) 

00  1112 

MRITE(Stlll) CA(J,I>fI»ltH,l) 

111  FORMAT  C"  **•  10 (  2X  t  E1Q«  4) ) 

1112  CONTINUE 
RETURN 
END 
C 

c 

C 

C 

SUBROUTINE  KPOWP CM,NP ,N, IOIH,RI 
C  FINOS  R»M**NP 

DIMENSION  M (IOIM, IOIM)  ,MI0IN,I0IN> 

DIMENSION  Rl(lf«lS)fR2(lt*lt) 

00  193  J-1.I0IM.1 
00  194  I*l» XDINt 1 
RlfI»J>*«.« 

R2(Xf J)*9«9 
194  CONTINUE 

Rif J»J)«1.I 
193  CONTINUE 

00  199  JJ»1,MP,1 

CALL  NULT (N»RltR2»N«N«Nv XOXNI 

CALL  COPY  (R2»Rl»  N«Ir<r  H) 

199  CONTINUE 

CALL  COPY (R2,R,N tIOIM) 

RETURN 

ENO 

C 

c 

C 

c 

SUBROUTINE  NOIIE  (RM  SNOTS,  OUTNEAN.MM) 

C  SUflTOUTINE  NOISE  CALCULATES  THE  VALUES  OP  THE  MEASUREMENT  NOISE 
C  COMPONENTS  USXMC  A  RAMON  NUMBER  GENERATOR  MODELLED  AS  SAUSSXAN 

CAUSS-f, 

00  331  X«l»12»l 
CAUSS»6A9SS*RAMF I Z1B) 

933  CONTINUE 

OAUSS« 6 AUSS-* . ♦OUTNEA N 
NM*GAUSt*RM SNOTS 
RETURN 
OH 
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SUBROUTINE  HGRAPHCX ,Y ,N, IC,NO,NP,NS) 

C  IF  I D(  COC  BOX  IN  UPPER  RIGHT  CORNER 

C  IS  NOT  PLOTTED 

01 HENSION  X(i),Y(l)  ,ID<1)  f  XF(N0.EQ.2)  GO  TO  30 
IF  CNO.LT.O)  GO  TO  1C 

CALL  SCALE(X,7.,N,i)  S  CALL  SC*LE(  Y,'  .  ,N,i) 

10  CALL  PLOT  (5.5,0.  ,-3)  1  CALL  PLOT  (C  .  ,11.  ,3) 

CALL  PLOT  (-1.35,1.35,  3) 

CALL- PLOT (-7.13, 1,35, 2)  S  CALL  PLOT (-7. 16 ,9.65 , 2) 
IF(I0<1)  .EO.CaO)  GO  TO  25 

CALL  PlOTC-7. 65,9.55,  3)  i  CALL  PLOT  C  -7.  C5,7.55,  2» 
00  29  1*1, 7, 2 

20  CALL  SYMBOL  (I*. 1-6. 9, 7. 85  , . 07,10 (I) ,90 . , 26) 

CALL  PLOT  (-7.05,7.5?,  3)  *  CALL  PLOT  (  -ft.  05,7. 55, 2> 

CALL  PLOT  (-6.05,  9.65,  2)  $  CALL  PLOT(-7.05,9.i>5,2» 

CALL  PLOT (-7.15,9.65, 3) 

25  CALL  PLOT (-1.35, 9.65, 2)  S  CALL  PL0T(-1.36,1.35,2) 
CALL  SYMBOL  (-6.6r-  ,1.15  ,.  1,10(13)  ,3  •«-*)) 

CALL  AXIS(-1«85,2.1*ID(9)»-2J*7.,93.  ,X  ( N*i )  ,  X  (N  *2)  ) 
CALL  AXIS  (-1.63, 2.1, 10(11),  20,5.  ,lSt.,Y(N*l>  ,Y(N42>> 
3d  Y(N*2)*-T(N42> 

X(N4l)*X(N41)-2,l*X(N42>  S  Y(N4ll*Y(N4l)M.85*Y(N42» 

CALL  LXNECY,X,N,l,NP,hS) 

X(N*1)  *X  (N41)  ♦£,  l*X  (N42)  %  Y(N41)»Y(N*1)-1.85*Y(N*2» 
Y(N42I*-Y(N42) 

RETURN  S  END 

SUBROUTINE  VGRAPH(X,Y ,N, ID,NO,NP,NS) 

OX HENSION  X(1),Y(1)  ,16(1)  $  IP(NO. EH.21  GC  TO  30 
IF  (NO.LT.O)  GO  TO  16 

CALL  StALF(X,«.9,N,l>  1  CALL  SCALE  (Y,7.C,N,  1) 

10  CALL  PL0T(8.5,3.,-3)  1  CALL  PLOT  CO.,  11., 3) 

CALL  PLOT  (-1.15,  1.35,  3) 

CALL  PLOT  (-7.15,  1.35,  2)  »  CALL  PL0T(-*M5 ,9.65,2) 

CALL  P;0TC-i.3J,9.5E, 2)  1  IF (10(1) «EQ. 90T)  GOTO  25 
CALL  PLOT (-1.A J, 9.»5, 3)  l  CALL  PL0T(-3. AS, 9.55, 2) 

00  29  X*  1,7,2 

20  CALL  SYNPOL  (-3. If ,5 .A -I* .16,. 07, 10 (I), 9. ,20) 

CALL  PLOT (-3. A f, 9. §5, 3)  i  CALL  FLOT( -3.45, 0 .55, 2) 

CALL  PLOT (-1.69, 8.3?, 2)  1  CALL  PLOT( -1.A5, 9.55,2) 

CALL  PLOT  (-1,35,  9,66,3) 

Cf  CALL  PLOT (-1.35,1.35, 2) 

CALL  SYMBOL (-6.0! ,1,19,. 1 ,10(13) ,6 ,,69) 

CALL  AXISC-t.*,  1.00  ,16(6), -20,6. 1,  1.  ,X(N4J  ) ,  XCN  +  t)) 
CALL  AXrS(-6.**, 1.85,10(11), 20, 7.1, 90  .,  Y  (N4i),Y(H42)» 
M  X(tt4lMXC«4i)46.*«Y(rt4t)  S  Y (M4i)  »V  (H41)  -1,05*Y  (N42) 
CALL  L?Nf(X,YvN,l,NPtNS> 

X(N4l)  at  (MU)  *6, -*X  (N4» 

HE  TURN  S  EMO 


I  YCN41)  *Y(N4l)4l,85»Y(H4{) 
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20.  ABSTRACT  (Continue  on  rover ee  eide  II  neceaeery  end  Identity  by  block  number) 

A  new  discrete  control  law  is  implemented  and  examined.  Robustness  to 
noise  and  model  errors  of  the  control  law  as  sample  rate  varies  is  analyzed. 
This  analysis  is  conducted  while  controlling  several  different  very  lightly 
damped,  single  Input/single  output  systems  which  are  representative  of  the 
flutter  dynamics  of  the  B-52E  wing. 

The  control  law  performance  1$  different  at  separate  sample  rates.  A 
distinct  range  of  sample  rates  are  found  to  have  a  better  response  to  noise 
than  other  sample  rates.  Another  range  is  found  to  be  more  robust  when  then 
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exists  an  error  in  the  models  used  to  calculate  the  closed  loop  control  law 
When  these  ranges  of  sample  rates  intersect,  the  robustness  characteristics 
at  those  sample  rates  is  found  to  be  good  with  respect  to  both  noise  and 
model  mismatch. 

As  theoretically  predicted  a  relationship  between  condition  number 
of  the  system  Hankel  matrix  and  robustness  seems  to  exist.  Hence,  these 
simulated  results  appear  to  validate  the  theoretical  results  on  robustness 
predicted  by  Reid,  but  on  the  other  hand,  these  simulated  results  indicate 
that  the  total  analysis  of  "robustness"  is  a  very  conplex  issue  and 
cannot,  at  this  point,  be  totally  predicted  by  such  a  parameter  as  simple 
as  the  cgndition  number  of  the  Hankel  matrix. 
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